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ABSTRACT 



Noise generation caused by quantization of signals in 
digital filters is considered using matrix state-space tech- 
niques. The results are expressed in terms of the correla- 
tion of the errors so that the assumptions made about the 
statistical properties of the noise sources can be evaluated. 
The effect of correlation, which is shown to be significant, 
depends upon the specific structure of the filter and the 
relative location of the error source both topologically and 
temporally . 

Specifically, the' results are applied to a newly extended 
set of canonic realizations of second order filters which are 
presented in the form of canonic arrays of coefficients. 

Output noise formulae for these realizations are presented 
for the two extreme assumptions of uncorrelated error sources 
and constant noise generation and are confirmed by experi- 
ment. These formulae are directly useful in analysis, design, 
and comparison. 

Additionally, new algorithms are presented, one for 
performing circular convolutions and the other for calculating 



the Discrete Fourier Transform. 
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I. INTRODUCTION 



A. IMPORTANCE OF DIGITAL SIGNAL PROCESSING 

One of the most pervasive advances in technology today 
has been in the area of digital signal processing. Since 
the time when the transistor made the proliferation of general 
purpose digital computers economically feasible, other ad- 
vances, particularly in integrated circuit technology, made 
practical the development of dedicated digital processors. 

From the common telephone and television, to rapid transit 
systems, manufacturing processes, speech therapy, and the 
exploration of space, digital signal processors have had 
their impact. In the Navy, digital signal processors are 
becoming more and more widely used in radar; sonar; missile 
and torpedo guidance; launcher control systems; combat 
information collection, dissemination and display; and of 
course communications systems. 

B. AREA OF INVESTIGATION 

A digital filter is defined [l] as a computational pro- 
cess or algorithm by which a digital (discrete-time and 
discrete-value) signal or sequence of numbers termed the 
input is transformed into an output digital signal. A 
digital filter may be implemented as a subroutine in a 
general purpose computer or as hardware In the form of a 
special purpose digital processor. 



Digital filters are a subclass of the set of discrete- 
time filters which also includes "sampled data" filters in 
which the signals may take on a continuum of values. The 
mathematical treatment of the discrete-time nature of these 
filters is common to both digital and sampled data filters. 
But in digital filters the discrete value property causes 
deleterious effects in the processing of signals. 

The reason that digital filters must operate on discrete 
valued signals is that there is only a finite number of 
digits available with which to represent the signal values. 
The only way to transform one set of numbers into another 
set of numbers is by a combination of the operations of 
multiplication and addition and it is also necessary in a 
digital filter that multiplier coefficients be representable 
in a finite number of digits. In the process of performing 
a multiplication, the length of the mantissa in the represen- 
tation of the product can be as large as the sum of the 
lengths of the mantissas of the multiplier and multiplicand. 
But in general that means the product is not representable 
by the number of digits available. Therefore, some method 
of quantizing must be performed during each sample time in 
order to maintain the representability of the signal quan- 
tities. This causes non-linear effects even in a filter 
which, when viev/ed only as a discrete-time filter, may 
appear to be linear in nature . 

Many authors have claimed that the errors generated by 
the quantization performed in the filter are uncorrelated 
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or at most negligibly correlated from sample-time to sample- 
time and between the sources of error. But when the phe- 
nomenon of limit cycles is considered it becomes obvious 
[Appendix A] that there is a great deal of correlation 
involved. Since the amount of correlation depends on the 
value of the coefficients [Appendix B] and the method of 
quantization, it follows that if there are more than one 
structure (algorithm) by which a particular filter may be 
realized and the values of the coefficients and the method 
of quantization are different between these structures, then 
the errors generated and their correlation will be different. 

The objective of this work is two-fold: (1) To formulate 

a general approach for investigating noise generation and 
propagation to include correlation and structural effects; 
and (2) to apply this approach to specific examples, partic- 
ularly the class of canonic second order filter structures 
which are derived herein. 

C. PREVIEW OP RESULTS 

After a summary of discrete time signal processing theory 
in Chapter II, an exhaustive search is made in Chapter III 
for all possible realizations of a second order digital 
filter in state variable form which are considered best in 
a certain sense. This search culminates in the conclusion 
that there are 58 ways to realise a given transfer function 
which meet the specified criteria. Thirty-eight of these 
realizations have not appeared previously in the literature. 
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In Chapter IV, following a method presented by Parker 
and Yakowitz [2,3], hut without assuming uncorrelated errors, 
a general expression for the output error variance for any 
order filter is derived. This expression is then applied 
to the second order case and then to the realizations pre- 
sented in Chapter III. When the assumptions of uncorrelated 
and highly correlated noise are made, the expressions for 
the output error variance reduce to two sets of formulae 
for these extremal cases. Examination of these formulae 
reveals for the first time the effects of structure on noise 
generation due to quantization. Appendix D verifies that 
these formulae yield an accurate expression for the output 
error properties. 

The quantization noise problem is particularly serious 
in recursive digital filters wherein the algorithm uses the 
results of previous calculations to generate present signal 
quantities. The fact that quantization errors are also fed 
back is the root cause of effects such as limit cycle genera- 
tion. On the other hand, non-recursive filters do not exhibit 
this effect; nor do recursive filters which produce no quan- 
tization errors which are fed back. In examining the nature 
of non-recursive and recursive finite impulse response (FIR) 
discrete-time filters in Chapter V, some interesting proper- 
ties of the signal quantities and algorithms are observed. 
These observations lead to an algorithm for computing con- 
volutions which is most efficient (in terms of number of 
operations performed) when completely overlapping blocks of 
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data from a continual data stream need to be processed in 
this way. Additionally a new method for computing the Dis- 
crete Fourier Transform (DFT) is presented which has advan- 
tages over other methods when a continual data stream is to 
be processed. These algorithms are performed without 
recursive error generation. 
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II. DIGITAL SIGNAL PROCESSING 



A. DISCRETE-TIME PROCESSING 1 

The most useful and most used tool in discrete-time 

signal processing is the Z-Transform. Even the Discrete 

Fourier Transform is a special case of the Z-Transform. 

This section reviews the applications of the Z-Transform, 

while the next section deals with the problems which arise 

2 

when this theory is applied to digital signal processing. 

1 . The Z-Transform 
a. Definition 

00 

Given a sequence (x(n)} the Z-Transform of 

n=- co 

the sequence is defined as 

00 

X(z) = Z[x(n) ] = Z x(n)z" n (2.1) 

n=~ co 

This definition is referred to as the two-sided Z-Transform. 
Another form of the Z-Transform is defined when x(n) = 0 
for n < 0. Then (2.1) becomes 

CO 

X(z) - Z x(n)z _n (2.2) 

n=0 



The use of the phrase discrete-time rather than digital 
is to emphasize that the theory in this section applies to 
"sampled data" signal processing, in which the data can take 
on a continuum of values, as well as digital signal process- 
ing, in which the data is restricted to discrete values [l]. 

p 

The object is not to reiterate known theory but rather 
to establish notation and definitions for use in later sections. 
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This definition is called the one-sided Z-Transform. Since 
(2.1) is more general, the term Z-Transform will refer to 
(2.1) unless otherwise specified. 

Both (2.1) and (2.2) can be viewed as Laurent 
series in the complex variable z. It can be shown [4] that 
for physical systems of interest (2.1) and (2.2) converge 
on the unit circle in the z-plane. Therefore, since the 
coefficients of a Laurent series are given by the integral 

x(n) = *C X(Z)Z " ¥ (2 ' 3) 



where C is a closed contour in the region of convergence 
and encircling the origin, (2.3) is taken to be the inverse 
Z-Transform and C is taken to be the unit circle. 

b. Relation to the Laplace-Fourier Transform 

When the sequence {x(n)> is the result of the 
equally spaced sampling of a piecewise continuous signal 
then the Z-Transform, X(z) is related to the two-sided 

A /S 

Laplace Transform, X(s), of the function, x(t), by 



X(z) 



z=e 



sT 



1 

T 



Z X(s 

k = -oo 



+ 



j 27Tk >. 
T ' 



(2.iJ) 



where j = "Y -1 and T is the spacing between samples. Note 
that this relationship is not one-to-one since there are an 

A 

infinite number of functions x(t) which can yield the same 

set of samples x(n). Mathematically, the mapping of the 

sT 

s-plane to the z-plane by z = e is many-to-one . Since 
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-1 -sT —1 

z = e , z is called the unit delay operator by the 
shifting theorem of Laplace Transform theory. 

Since the Fourier Transform can be equivalently 
viewed as the two-sided Laplace Transform evaluated at s = jw 
then the Z-Transform evaluated on the unit circle is related 
to the Fourier Transform by 



X(z) 



z=e 



j 



u>T 



i I X(j((0+^))=i Z X(J(w+ku> s )) 



where oj is the sampling frequency in radians. Notice that 

S 

HCe^ 1 ^) is periodic in w with period . 

2 . Functional Transforms 

In continuous linear theory, the relationship between 

A A 

the input, u(t), and the output, v(t), is given by 



V(s) = H(s)U(s) (2.6a) 



and 



v ( t ) 



/ h(x)u(t-T)dx 

— 00 



✓s A 

/ h(t-x)u(x)dx (2.6b) 

_oo 



A A A 

where V(s) and U(s) are the Laplace Transforms of v(t) and 

A A A 

u(t) respectively, H(s) is the Laplace Transform of h(t) . 

A 

The function h(t) is called the impulse response of the 
filter or equivalently the Green's function of the linear 

A A 

integro-dif ferential operator which relates v(t) to u(t). 
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Equation (2.6b) is the convolution of h(t) with u(t) denoted 

^ ^ A 

by h(t)*u(t). H(s) is called the transfer function or fre- 
quency response of the filter. 

In discrete-time linear filter theory, the 
analogous equations hold, i.e. 

V(z) = H(z)U(z) (2.7a) 

00 oo 

v(n) = E h(m)u(n-m) = E h(n-m)u(m) (2.7b) 

m=-» m=- 00 

where V(z), H(z) and U(z) are the Z-Transforms of v(n), 
h(n) and u(n) respectively. The sequence(h(n) } is called 
the impulse or unit pulse response and H(z) is the transfer 
function of the filter. 

A /\ 

To be physically realizable h(t) and h(n) must be 
zero for t < 0 and n < 0 respectively. This realizability 
condition is necessary in order that the filter not react 
to something which has not yet occurred. Thus this condition 
is referred to as causality. 

The transfer function, H(z), can always be written 

A 

as a ratio of polynomials in z, and H(s) as a ratio of poly- 
nomials in s . 

In designing a discrete-time filter, it is often 
the case that one tries to duplicate a well-known continuous- 
time filter. It is then possible to transform the transfer 
function in the s-domain into a transfer function in the 
z-domain. Several functional transform methods have been 
used [5]. 
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a. Standard (Impulse Invariant) Z-Transform 

The Standard Z-Transform maps the transfer func- 
tion from the s-domain to the z-domain by the transformation 
s T 

z = e the value of T and thus to having previously been 

s 

decided. This method always results in the unit pulse re- 
sponse of the discrete-time filter being identical to the 
samples of the continuous-time filter taken at intervals 
T of the independent variable, t, starting at t = 0, hence 
the equivalent name. Impulse Invariant Z-Transform. There 
is a difficulty in this method, however. By Equation (2.5) 
it can be seen that the frequency response of the discrete 

time filter for -to /2 < to, < co /2 will involve the sum of 

s — d — s 

the values of the frequency response of the continuous-time 

filter at to = to + kto , _« < ^ < “. This effect is called 

d s 

i coT ^ 

aliasing and will result in H(e J ) not being equal to H(jto) 
in the region of interest. They may be approximately equal 

A 

if H(jto) is very small outside half the sampling frequency, 
but no time limited function (which an impulse response is) 

A 

can be limited in frequency . On the other hand, when H(jto) 
is not small outside half the sampling frequency, H(e^ t0 ^') 
will bear no resemblance to the desired frequency response . 

For this reason other functional transforms are needed. 

b . The Bilinear Z-Transform 

Another method of transforming a transfer function 
from the s-domain to the z-domain is by the mapping 
z = (1 + ~)/(l - ~) whi ch has the inverse mapping 
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2 z — 1 

s = T^ z + 1 ~) * mapping, however, distorts the frequency 

response. Therefore it is necessary to counter-warp the 
desired radian frequency response before applying the trans- 
formation. This still does not yield an exact equivalence 
between the two frequency responses because the pole-zero 
locations are not matched. 

c. The Matched Z-Transform 

The Matched Z-Transform is also a Bilinear 
Z-Transform using the same bilinear mapping as before. But 
the way the transfer function is modified is to replace each 

A 

pole or zero of H(s) given by 



s 




0 



(2.8a) 



with a pole or zero of H(z) given by 



s . T 

z-z^ = z- e 1 



(2.8b) 



This process effectively prewarps both the radian frequency 
response and the neper frequency response, 
d. Higher Order Transforms 

Examination of the effect of the preceding func- 
tional transform techniques reveals that the algorithm implied 
by the transfer function in the z-domain is a numerical inte- 
gration routine where, for the Standard Z-Transform, the 
routine is of zero— order (rectangular) and for the Bilinear 
Z-Transform it is first-order (trapezoidal). It is possible 
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then to hypothesize using higher order numerical integration 
algorithms for the transformation between the s- and z- 
domains . This would have the effect of raising the order 
of the transfer function (which is the order of the denom- 
inator polynomial) in the z-domain. This increases the 
complexity and cost of the filter and thus is not usually 
done explicitly in practice, but may be done implicitly as 
observed in the next section. 

e. Finite Impulse Response (FIR) Filters 

When a filter’s impulse response (which is zero 
for n < 0) is also zero for all samples after some sample 
h(N), the filter is called a finite duration impulse response 
or simply finite impulse response (FIR) filter. In this 

. . M 

case, the denominator polynomial of H(z) is z" and the 

numerator is I h(n)z N_n . If an FIR filter is being used 
n=0 

to approximate a continuous-time filter of order less than 
N, then some combination of higher order integration routines 
are implicitly being performed. But the intent in designing 
an FIR filter is usually to approximate some property of 
the desired frequency response such as linear phase [6], 
minimum phase [73, or magnitude response [8,9] often without 
reference to the pole-zero locations of a continuous-time 
filter and letting the order of the filter be that which is 
v needed to meet some design criterion such as minimax or 
mean-square difference. 
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3. Transfer Function Realizations 



Once the functional form of the filter has been 
decided upon, it is then necessary to design the realization 
of the algorithm to be performed, 
a. Definitive Forms 

The most direct way to realize a transfer function 
is by performing the operations implied by that function. 

(This is called the direct realization for obvious reasons.) 
The transfer function is usually written as the ratio of 
two polynomials in z, such as 



K ( z ) 



M 

Z a. z 
i=0 x 
N 

Z b . 
0=0 0 



(2.9) 



with N _> M for realizability." 5 Then it is possible to 

-N 

multiply numerator and denominator by z to yield 



H ( z ) 



M 

E a . z 
i=0 1 
N 

Z b .z 
j=0 J 



i-N 



0-N 



N 
Z 

i=N-M 
N 

Z b„ - z 

o’=o N - J 



a N-i Z 



-0 



-l 



V(z) 

uTzT 



( 2 . 10 ) 



o 

It will be shown in Chapter III, that in some contexts 
it is necessary that N be greater than M for realizability. 



26 



Recalling that z ^ is the unit delay operator, 

■“k 

the product z X(z) is the Z-Transform of x(n-k) . Therefore, 
by cross multiplying the last equality in (2.10), taking 
the inverse Z-Transform of both sides term by term and 
letting b^ be unity (without loss of generality), the time 
domain equivalent of (2.10) is given by 

N N 

v(n) = 2 a.. . u(n-i) - Z b„ ,v(n-j) (2.11) 

i=N-M j=l J 

Equation (2.11) is realized directly in Figure (2-1) which 
has 2N unit delays indicated by z”\ But by defining an 
intermediate variable x(n), such that 



N 

x(n) = u(n) - Z b N _jX(n-j) 



it can be shovm that (2.11) becomes 



( 2 . 12 ) 



N 

v(n) = Z a. T . x(n-i) (2.13) 

i=N-M rg " 

Equations (2.12) and (2.13) represent the algorithm of 
Figure (2-2) which is known as a canonic realization. It 
is canonic in that it has the least possible number of delay 
elements (N) . 

b. Reduction to Lower Order Forms 

VJhen the numerator and denominator polynomials 
are known in factored form, it is possible to realize the 
transfer function as a combination of lower order transfer 



functions . 
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(1) Parallel Realization by Partial Fraction 
Expansion. If the transfer function has simple poles, then 
it can be written as the sum of first and second order trans- 
fer functions (sometimes called sections) by performing a 
partial fraction expansion. These can then be realized 
individually and then placed in parallel between the input 
and output as in Figure 2-3. If the transfer function has 
multiple poles, higher order sections will be required and 
will generally result in a higher number of delays and 
multipliers than in the canonic realization. 

(2) Cascade Realization by Factorization. In 
the factored form there is a definitive way to realize the 
overall transfer function by arbitrarily associating zeroes 
and poles and writing the transfer function as a product of 
lower (usually first and second) order sections as in Figure 
2-H. 

(3) Hybrid Realization. The operations of (l) 
and (2) above may also be used in many different combinations 
to form many different hybrid realizations, such as in 
Figures 2-5a and 2-5b. 

c. Canonic First and Second Order Sections 

Since any transfer function can be realized as 
a combination of first and second order sections, these low 
order forms take on a great deal of importance in the design 
of discrete time filters. Figures (2-6a and b) show the 
canonic realizations of first and second order filters. 
Jackson [10] discussed variations of Figure (2-6b) which are 



28 



just as efficient in terms of operations. Hess [11] 
extended these to 24 forms. This work extends these to 
58 forms. 



d. Other Forms of Realizations 



There are several other forms for the realization 



of a transfer functions . Among these are the so-called ladder 
structures [12,13] found by continued fraction expansions 
and the Fettw r eis structures [l4] which are found by direct 
analogy to RLC ladder structures . Both of these types are 
ladder structures so care in terminology is necessary to 
avoid confusion. 

4 . The Discrete Fourier Transform 

As mentioned in the first paragraph of this section 
the Discrete Fourier Transform (DPT) is a special case of 
the Z-Transform. The DFT is the Z-Transform of a sequence 
of length N evaluated at N equally spaced intervals around 
the unit circle in the z-plane and can thus be defined by 



X(kfl) = X(z) 




0 < k < N-l 



(2.14) 



where 



0 = 2ir/NT 
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The Inverse Discrete Fourier Transform (IDFT) is given by 

, N -i , 

x(n) = rr I X(kfi)W nK 0 < n < N-l (2.15) 

1N k=0 

Notice that if n and k are outside the limits [0,N-1] then 
(2.1*0 and (2.15) are periodic with period N. 

The IDFT of the product of the DFT’s of two sequences 
of length N is by definition reducible to the circular con- 
volution of the two sequences. The process of circularly 

2 

convolving two sequences of length N requires N multipli- 
cations, while the use of the DFT to accomplish the same 

2 

operation takes 3N + N multiplications so it is more effi- 
cient to perform the convolution directly. There were many 
early attempts to make the DFT more efficient, including 
the Goertzel algorithm [4], Then Cooley and Tukey [15] 
found a way to reduce the number of operations in finding 
the DFT which has become known as the Fast Fourier Transform 
(FFT) [16-19]. 

There are many variations of the algorithm now, 
including algorithms for N a composite number of mixed radix, 
algorithms which require sorting in time (decimation in time) 
as in the Cooley-Tukey algorithm and algorithms which require 
sorting in frequency now known generically as Sande-Tukey 
algorithms . Further reductions in the number of operations 
required has been accomplished by observing symmetries and 
common factors. 
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The DFT has also been widely used in the discrete 
filter design process [8, 20]. 

5 . Other Transform Techniques 

Other transform techniques utilized in discrete time 
signal processing include the Chirp Z-Transform [21], Fermat 
Number Transforms [22, 23], and the Walsh-Hadamard Transform 
[2^,25] and the Hilbert Transform [26], 

B. DIGITAL PROCESSING 

When the theory of Section A is applied to digital 
processing, the discrete valued nature of the numerical 
representation required produces nonlinear effects which 
disturb the solutions to what otherwise would be linear 
equations. Among the problem areas of concern are analog 
to digital (A/D) conversion errors, coefficient representa- 
tion (the approximation problem) and quantization errors 
involved in arithmetic operations . 

1 . A/D Conversion 

The problem of A/D conversion is not one which lends 
itself to solution. It is a problem which must be contended 
with, the only way of reducing it being to use more signifi- 
cant digits in the digital signal representation thus making 
this purely a cost vs. error tradeoff consideration. Bennett 
[27] has presented spectral studies of A/D conversion noise, 
thus allowing the designer to represent this error as an 
injected noise at the input of the digital processor. Often 
the quantization level of A/D converters is much larger than 
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that of the digital processor [28], The processor may then 
include some filtering action to reduce the effect of the 
A/D conversion errors . 

2 . The Approximation Problem 

When attempting to synthesize a digital filter as 
an approximation of some required characteristic such as 
impulse response or frequency response, the determination of 
filter coefficients is complicated by the finite word length 
representation allowed. This causes error minimization 
procedures to become iterative nonlinear routines. 

The sensitivity of digital filter algorithms to 
shifts in the ideal coefficient values has been studied 
extensively. Kaiser [29] concludes that for a direct filter 
realization, the sensitivity of pole positions increases with 
the order of the filter. This has been confirmed by Knowles 
and Olcayto [ 30 ]. Rader and Gold [31], in studying the 
sensitivity of second order filters, conclude that two 
coupled first order sections are less sensitive than a single 
second order section. But it should be observed that two 
coupled first order sections require more multipliers than 
canonic forms (see Chapter III) and are merely a general 
state variable realization which is not canonic. Mantey 
[ 32 ] studied state variable representations and concluded 
that digital filters should be realized as parallel or cas- 
cade combinations of first and second order sections . This 
is consistent with the conclusions of Ka3 ser above . Ladder 
structures in general have been found to be less sensitive 
to coefficient inaccuracy [33]. 
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3 . Quantization Errors 



Arithmetic operations in digital filters almost 
Invariably lead to the requirement to requantlze the result 
of the operation at some time during one cycle time. Mul- 
tiplication increases the word length of the signal. Addi- 
tion and multiplication may exceed the dynamic range of the 
numerical representation. Floating point addition may re- 
quire quantization before the operation since both the sum- 
mand and addend must have the same exponent. 

The choice of arithmetic mode significantly effects 
the level of quantization noise. Since fixed-point arith- 
metic is easier to implement it is most popular. But the 
relative noise level can be severe if the full dynamic range 
of the registers is not being used effectively. 'This requires 
scaling the driving function. Floating point arithmetic 
utilizes the exponent for automatic scaling. An intermediate 
method is known as block floating point [3*0 where the expon- 
ent for a whole block of data (such as a single second order 
section in a cascade realization) is maintained in a single 
register . 

Oppenheim and Weinstein published a review [351 °f 
the effects of quantization errors in digital filter and 
Fourier Transform algorithms. This review includes the work 
of Liu and Kaneko [36] on floating point arithmetic, Knowles 
and Edwards [37] and Gold and Rader [38] on fixed point 
arithmetic, and Weinstein and Oppenheim [39] where all three 
modes were compared on the basis of signal-to-noise ratio. 
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It was determined there that, in general, floating point 
arithmetic was best with block floating point falling 
between the other two modes. 

The choice of quantization mode also effects quanti- 
zation. For addition, two's complement arithmetic exhibits 
disastrous errors when overflow (upper level quantization) 
occurs [40, ill]. Zeroing and saturation arithmetic yields 
better results [42], For multiplication and addition, 
lower level quantization methods include rounding, trunca- 
tion and sign-magnitude truncation. Jackson [43] has studied 
the interaction of upper level and lower level quantization 
errors . 

Since quantization error is a nonlinear effect, 
experience has shown that the phenomenon of limit cycles 
exists in digital systems just as in continuous ones. Limit 
cycles have been extensively studied both in statistical 
and deterministic ways . 

Certain bounds on the magnitudes of limit cycles 
in particular and accumulated roundoff errors in general in 
second order sections have been proposed by many authors and 
are tabulated in Table 2-1. 

These bounds, however, have been derived assuming 
uncorrelated error sources, a restriction which will not 
be assumed in the derivations of Chapter IV. 
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Author 



Reference 



Bound 



Condition 



Jackson (estimate) 44 




Jackson (bound) 10 




Bonzanigo 45 


■ r^i f“0 f /2 

1- af+b u > I s / ^ 


Sandberg and Kaiser 46 


2 n/2 b>0 

( 1-b ) [ 1- (a / 4b ) ] 1/2 

i i to 

l a l— 1+b 

0<f / f s <| 


I-JaT+b' b -° 

1 1 to 

l a l- 1+b 


Long and Trick 47 


1- a + b 


1 - *b 2 ( 2 -i)<a 2 <i.b 

(l-Vb) Vb 

b>0 


1+ -/b 2 ,.,2/ 2 , a 

: = 75 - a <4b (-=-1) 

(l-b)[l-(aV4b)] V ^ Vb 

b>0 


Yakowitz and Parker 2 


l-l'al+b /(| 


J 7 J b>a 2 /H 

(1-b ) [1- (a /4b )] /d 



f =saupling frequency ,1+az ^+bz ‘"-denominator of H(z) 

s 

TABLE 2-1 Present Bounds on Quantization Error 
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In statistical studies of quantization errors, the 
expressions for the output error variance, assuming uncor- 
related errors and zero mean, uniformly distributed statistics, 
usually involve contour integrations of the form [4,48] 



= Z 
i 



2irj 



15 G i (z > G i<?> f 



( 2 . 16 ) 



2 

where Oq is the output error variance 

X. T_ 

and G^(z) is the transfer function between the i v error 
source and the output. This form of analysis requires the 
determination of all the G^(z) and the performance of all 
the contour integrations. Furthermore, the assumption of 
uncorrelated error sources provides no insight into the 
effects of structure on the output error variance. Nor does 
this method allow the effect of assumptions made about the 
distribution and correlation of the errors to be examined. 
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p 

H(z) = K + S G.(z) 
P i=i 1 



Figure 2-3. Parallel Realization of H(z) 




H(z) = K n G.(z) 
i=l 1 



Figure 2-4. Cascade Realization of H(z) 
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H(z) = K^Kg + K 3 G 1 ( 2 )G 2 (z) + K 4 G 3 (z)G 4 (z)] 



Figure 2-5a. Hybrid Realization cf H(z) 




H(z) = K n [K 2 + H -j ( z ) + H 2 (z)][K 3 + H 3 (z) + H 4 (z)] 



Figure 2-5b. Hybrid Realization of H(z) 
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Figure 2-6a. Canonic First Order Section 
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Figure 2-6b. Canonic Second Order Section 
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III. CANONIC REALIZATIONS OF SECOND ORDER SECTIONS 



A. INTRODUCTION 

In section IIA3b it was shown that any rational transfer 
function could be expressed as a parallel, cascade or hybrid 
parallel-cascade combination of first- and second-order 
transfer functions ( or -sections ) . For reasons of noise 
generation, coefficient sensitivity, economics, logistics 
and others, these low order sections assume great importance 
in the field of digital filtering. 

Two ways of realizing a second order section (by the 
direct or the so-called canonic realization (see section IIA3c) 
have already been given. The "canonic" form has been and, 
in general, still is the way second order sections are 
realized in practice. However, when one considers the second 
order section in its state variable form, one realizes that 
there are an infinite number of solutions depending upon 
the manner in which the basis is represented in the state 
space. And, in fact, many of these are canonic in the sense 
of minimizing the number of operations required. This mini- 
mization also implies a minimum number of noise error sources 
due to quantization after multiplication since the number of 
multiplication operations is minimized. 

In this chapter, some inquiries are made into the relation- 
ship between the transfer function and the state and output 
equations of a second order section. These relationships 



will then be exploited to find a group of state and output 
equation sets considered best in the sense of minimizing the 
number of quantization error sources. These sets are reduced 
to a number of "arrays" for compactness of notation. Half 
of the realizations represented by these arrays have not 
appeared in the literature to date nor has such an exhaustive 
search of allowable state equations been performed. The 
next chapter will delve into the relative noise properties 
of the realizations derived in this chapter with a view 
toward selecting optimal structure. 



B. . THE TRANSFER FUNCTION 
1 . Basic Forms 

The most general form of the second order transfer 
function is given by 



H ( z ) 



x -lx -2 
a Q + a^z + a 2 z 

b 0 + V 1 + b 2 z' 2 ’ 



V(z) 

uTzT 



( 3 . 1 ) 



where v(n) is the output and u(n) the input. 

In this form, it appears that there are six coeffi- 
cients and hence six sources of error due to quantization 
after multiplication. It is also possible to write equation 
(3*1) in the form 



H(z) 



-1 X -2 
ctg + ct-^z + a^z 

zr r? 

1 + az + bz 



( 3 . 2 ) 



f l3 



where 



a 0 ' ( 3 « 2a) 

“l = a l /b 0 (3.2b) 



a 2 ~ a 2 //b 0 



a - b x /b 0 



b = b 2 /b Q 



(3.2c) 
(3. 2d) 
( 3 . 2e ) 



so that there are only five non-unity coefficients. (Note 
that it is not necessary to require the condition b^ j* 0 
since if b Q we re zero the denominator of (3.1) would not be 
a second order polynomial in z 3 hence H(z) would not be a 
second order transfer function.) So H(z) has been written 
in (3*2) in a form that has at most five quantization error 
sources . 

From the form (3.2) it is also possible to rewrite 

H(z) as 



H(z) = K 



d + 



-1 




“2 


cz 


+ 


e z 


1 + az" 


•1 


+ bz ~ 2 j 



( 3 . 3 ) 



where 



“o * 0 



1 a 0 * 0 



0 a 0 = 0 



a o 



I lij 



d 



0 



0 



(3.3a) 
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(3-3b) 
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(3.3c) 
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a 0 ¥ 0 



a„ = 0 
0 



(3.3d) 



by dividing the denominator of (3.2) into the numerator in 
ascending powers of z~ 1 . 

Another form is possible if b ( or equivalently b 2 ) 
¥ 0. Then the transfer function may be written 



H(z) = K' 



d' + 



, , , -1 
c 1 + e 1 z 

1 + az" 1 + bz~ 21 



(3.4) 



where 



d' = 




a 2 ¥ 0 



a 2 = 0 



(3.4a) 
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t = 



“o b 

C(„ 



(X, 



- 1 



a 2 / 0 



a 2 = 0 






a 0 b 2 

a 2 b 0 



- 1 



u 0 



a 2 ¥■ 0 



a 2 = 0 



(3.4b) 
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aT 



a i 
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a 2 ^ 0 



a 2 = 0 
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a i b 2 

a 2 b 0 



a 2 / 0 



a 2 = 0 



(3.4c) 



cx. 



v t = 



a 2 ^0 



a 2 = 0 



a. 



= / 



a 2 ^ 0 



a 2 = 0 



(3.4d) 



by dividing the denominator of (3-2) into the numerator in 
descending powers of z ~^ . Note that the condition b ^ 0 
is equivalent to saying that there is no pole of H(z) at 
the origin in the z-plane . If there is^ a pole at the origin, 
then the other pole is real and two first-order sections are 
recommended rather than one second-rder section. So, in 
any practical second-order section it will be considered 
true that neither b Q nor b 2 is zero. 

It is still true, in the forms of (3-3) ana (3-4), 
that H(z) can be realized with at most five quantization 
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error sources. However, the coefficients K and K* need nor 
be considered part of the transfer function since they are 
scaling factors. Most often, in the practical realization 
of higher order transfer functions by cascade, parallel or 
hybrid forms (or even in realising a second order system) a 
scale constant will be used which is not K or K' but rather 
some number which, hopefully, will prevent overflow in the 
adders. Any gain factor necessary to restore the system to 
some overall gain factor is then included at the final out- 
put. Considering only the terms within the parentheses of 
equations (3*3) and (3.4), then there are at most four 
quantization errors. 

If it is the case that a^ = 0 (implying an infinite 
zero of H(z)) or = 0 (implying a zero at the origin), 
then it is possible to factor a scaling constant out of the 
rational part of (3-3) or (3.4) respectively, further reducing 
the number of error sources to at most three (four if the 
scale factor is included) . 

Or it could be that one or more of the coefficients 
is an integer in which case there can be even fewer error 
sources. But this case is not sufficiently general to 
warrant special consideration in developing algorithms 
for broad application or mass production. 

In the development that follows, four cases will be 
considered general enough for inclusion. These are the two 
cases of equation (3*3) with d = 0 or 1 and the cases of 
equation ’with d* = 0 or 1 . 
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It will be seen that equation (3.1) or (3-2) realized 
in the '•canonic" way in section IIA3c is identical to one of 
the canonic arrays to be developed, thus showing that it is 
not an additional realization. 

2 . Realizability 

There is a subtle fact regarding the physical reali- 
zability of a digital filter which is not readily apparent 
in the mathematical treatment of digital signal processing. ' 
Because the operations involved in performing an algorithm 
cannot be accomplished instantaneously, there must be an 
inherent delay between the time that the input arrives and 
the time the output is available for further processing. 

If the filter is being used at very low sampling rates com- 
pared to the operation speed of the devices, this delay may 
seem negligible but it exists nonetheless. On the other 
hand, if one is to design a filter that is usable in a broad 
field of application it is useful to consider the highest 
possible sampling rate. Then the magnitude of this delay 
is on the order of z - ^. 

When the difference equation associated with equation 
(3.2) is examined, the present output, 

v(n) = a^u(n) + a^u(n-l) + a 2 u(n-2) - av(n-l) - bv(n-2), 

(3.5) 

is seen to occur at the same time as the present input. But, 
by the previous discussion, this is impossible in a physical 
sense. Therefore we say that equation (3-2) is not, strictly 
speaking, realizable. 



From ( 3 - 5 ) it is seen that the factor which violates 
the condition of realizability is a^uCn) . This implies that 
(3.2) is not realizable because the coefficient, a^, in the 
numerator does not have a delay associated with it. Another 
way of saying this is that if the numerator and denominator 
of H(z) are multiplied by the lowest power of z which will 
clear H(z) of all factors, z”^, then the power of z multiplying 
c(q is equal to the highest power of z in the denominator. 
Therefore, in order to be realizable, any transfer function, 
H(z), written in terms of positive powers of z, must have 
a numerator polynomial of lower order than the denominator 
polynomial. If H(z) does not meet this criterion, it is 
always possible to realize z~ k H(z) where k is an integer 
large enough to force the transfer function into this condi- 
tion. In the case of equation (3.2), k = 1 will make the 
transfer function realizable when cXq ^ 0 . It may be argued 
that this effectively increases the order of the transfer 
function. This is of little consequence since the time 
representation of the output will be identical except for 
a shift of k sampling intervals. 

On the other hand it is not always necessary to make 
a second order section realizable in this sense . If the 
second order section is part of a cascade realization of a 
higher order transfer function, for example, it is only 
necessary to ensure that the overall transfer function is 



realizable . 



The question of realizability will be used in the 
analysis of possible sets of state equations in the next 



section in order to reduce the number of sets to be considered. 
C . STATE EQUATIONS 

The state equations for a second order digital filter 
consist of two first order difference equations involving 
the states and input (s) of the filter. The output equation(s) 
describes a combination of the states and input (s) which will 
be seen external to the filter as the output(s) . Since a 
transfer function is a relationship between one input and 
one output, only state and output equations with one input 
and one output will be considered. 

1. Basic Forms 



Define the vector, x_(n) , to be the vector of states 



in a second-order digital filter at time n, so that 




x(n) 



( 3 . 6 ) 




Now define the state and output equations of the 



filter as: 



x(n) = Ax(n-l) + Bu(?) 



(3.7a) 



v(n) = Cx(?) + 6u(?) 



(3.7b) 
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where 



A is a 2x2 matrix 
B is a 2x1 matrix 
C is a 1x2 matrix 
6 is a lxl matrix or scalar. 

The question marks in parentheses indicate an uncertainty, 
at this point, of the proper time index to be used. The 
n on the left side of equations (3*7) has been established 
as a reference point from which to determine the other time 
indices. Since it is already known that equation (3.7a) is 
a set of first order difference equations, it follows that 
A must multiply x(n-l) thus establishing that time index. 

Suppose the time index of u in (3.7a) is n-p, that 
of x in (3.7b) is n-q and of u in (3. 7b), n-r, where p, q, 
and r are integers greater than or equal to zero (since 
indices implying future data are not allowed) . Then 

x(n) = Ax(n-l) + Bu(n-p) (3.8a) 

v(n) = Cx(n-q) + 6u(n-r) (3.8b) 

Taking the Z-transform of equation (3.8a) and solving 
for X(z) yields 

X(z) = [I - Az -1 ] -1 BU ( z ) z“ P (3.9) 
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where I is the identity matirx. Taking the Z-transform of 
equation (3.8b) and substituting (3.9) for X(z) yields 

V(z) = CX(z)z“ q + 6u(z)z“ r 



transfer function between U(z) and V(z), which can be related 
to the transfer function in question. Two things need to 
be determined, (1) the values of p, q and r and (2) the 
entries in the matrices A, B, C and <5 . The first shall be 
taken up here; the second will be determined in the next 
section . 



= {C[I-Az 




(3.10) 



The factor inside the braces in (3.10) is thus a 



The quantity C[I-Az - ' 1 '] "'"B in (3.10) is a scalar 
formula in z ~ ^ given by 



$(z) = C[i-Az _1 ] 1 B 




Y + £2 



(3.1D 



where 




(3.11a) 




(3.11b) 



a = + a 22^ 
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^ a H a 22 " a i2 a 21 



(Note, b^ and b^ here are the entries of the matrix B, not 
the denominator coefficients of equation (3.1).) If equations 
(3.8) and thus equation (3.10) are to correspond to the trans- 
fer function H(z), then H(z) has to be of the form. 



H 



(z) = 6z _r + $(z)z~^ +c ^ 



- ,_-r + (Y+ez-bz-^ 
z T TT ^2 

1 + az + Bz 

. 6z -r + aSz- (1+r) + S«z- (2+r) + Y z- (p+q) + EZ -<P +q+1 > 



i x -r v . -2' 
1 + az + bz 



(3.12) 



If the requirement is now placed on the numerator 
of H(z) to be at most quadratic in nature (i.e., disregarding 
overall delay, the difference between the highest and lowest 
power of z ^ is no greater than two) it is possible to restrict 
the allowed values of p, q and r. 

First of all, the only indices available in the 
states of the filter are n and n-1, so q can only be 0 or 1. 

Nov; the highest power of z ^ in (3.12) is either 
p+q+1 or r+2 and the lowest is either r or p+q . The following 
logical arguments will place further conditions on p and r. 
a. argument (A) 

If p+q+1 >_ r+2 and p+q r then r-q+1 <_ p <r-q . 
But this is a contradiction, so this case can never arise. 
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b. argument B 

If p+q+1 >_ r+2 and r <_ p+q, then the stronger 
condition is that p >_ r-q+1, but the quadratic requirement. 
(p+q+1) - r <_ 2, implies that p £ r-q+1, so under these 
conditions p = r-q+1. Then for q = 0, p = r+1, and for 
q = 1, p = r. 

c . argument C 

If r+2 >_ p+q+1 and p+q £ r, then the stronger 
condition is that p ^ r-q, but the quadratic requirement, 

(r+2) - (p+q) £ 2, implies that p ^ r-q, so that p must be 
such that p - r-q. Then for q = 0, p = r, and for q = 1, 

P f r-1. 

d. argument D 

If r+2 _> p+q+1 and r <_ p+q then r-q <_ p <_ r-q+1 
and the quadratic requirement, (r+2) - r <_ 2 is automatically 
satisfied. Then for q = 0 , p = r or r+1, and for q = 1, 
p = r or r-1. 

Arguments B through D taken together imply that p 
and r differ by at most 1. Recalling that p and r are 
associated with both indices of u and only u in equations 
(3.8), it can be assumed without loss of generality that 
p and r take on only the values 0 and 1 since if they are 
greater than that, the signals u(n-p) and u(n-r) can be 
considered delayed versions of some other signal u'(n) or 
u'(n-l). Then the signal u'(n) or u'(n-l) could be used 
as the input . 
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Table 3-1 lists all possible combinations of p, q 
and r such that they take on only the values 0 and 1, along 
with the associated transfer function. In addition, the 
values of p+q, p+q+1 and r+2 are tabulated along with the 
appropriate logical argument which applies to their relative 
sizes. Referring to those arguments and the corresponding 
values of p, q and r, it is seen that lines (b) and (g) 
of Table 3-1 violate the conclusions of the argument with 
respect to the relative values of p and r. Thus the trans- 
fer functions of lines (b) and (g) cannot be related to a 
second order transfer function whose numerator is at most 
quadratic in nature. 

The state equations and output equation of line(h) 
in Table 3-1 can be seen to be identical to those of line 
(c) except that the input in both equations of (h) has been 
delayed one sample period. Therefore the internal structure 
would be the same except for a delay in the input line. 

The equations of line (e) can also be seen to be identical 
in nature to line (h) by rewriting the output equation of (e) 
as v(n-l) = Cx(n-l) + 6u(n-l) instead of 

v(n) = Cx(n) + 6u(n), whereas in (h) , v(n) = Cx(n-l) + 6u(n-l). 
In other words , the output of the equations for line (h) 
is a delayed version of the output corresponding to line (e), 
so the structure will be the same except for a delay in the 
output line . 

Similarly the equations of (f) represent a delay 
added on the input line of (a) and the equations of (d) 
represent a delay on the output line of (a) . 
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These two groups of essentially identical state and 
output equations are indicated as groups S and T respectively 
in Table 3-1. (Note that for 6=0 any of the transfer 
functions in Table 3-1 could belong to either group, including 
the two that apparently violated the logical arguments since 
for 6=0 the value of r is irrelevant. In order to remain 
as general as possible, 6 will remain arbitrary.) Thus the 
transfer functions associated with each member within a 
particular group represent essentially the same transfer 
function with the possible exception of an overall delay. 

For simplicity, two transfer functions are chosen to 
represent their respective group, namely 

Hg'(z) = [6 + *(z)z“ 1 ]z" 1 (3.13) 



and 



H t (z) = [6« + <Kz)]z 1 (3- 1*0 

where $'(z) has the same form as $(z). Equation (3.13) is 
taken from line (h) of Table 3-1 and (3.1*0 from line (f). 
The primes in (3.1*0 have been added for purposes which will 
become clear in the next section. The reason for the group 
designations and subscripts, S and T, is historical and will 
be made clear at the end of this section. These two forms 
have been chosen because they have an overall delay and are 
thus physically realizable . 
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It is now possible to make some specific connections 



between the two forms of H(z) given in equations (3.3) and 
(3«*0 and the state equations associated with (3.13) and 
(3-1*0 respectively. 

D. THE COEFFICIENTS OF THE STATE EQUATIONS 

Given a transfer function in the form of equation (3.3) 
(with an overall delay added for realizability) and a set of 
state and output equations which yield a transfer function 
in the form of (3.13) , the two systems can be made equal by 
setting their transfer functions equal and solving for the 
appropriate coefficients." 1 " In this case (disregarding 
scaling constants) 

cz + ez“^ -1 L , yz + ez~ 

IT Zp z = 6 IT 

1 + az + bz j l l + az + 3z 

= H s (z) (3.15) 

Equating coefficients, it is found that for (3-15) to be 
compatible 

a = a = -(&-)! + a 22^ (3.16a) 





^The derivation which follows is essentially due to Parker 
and Hess [A 9 ] which fellows a similar effort by Jackson [ 1 0 j - 
The derivation is reiterated here for continuity . The only 
difference is some minor changes in notation and the use of 
realizable forms of the transfer function. Parker and Hess 
used equation (3.3) without an overall delay and the state 
equations corresponding to line (e) of Table 3-1 with a change 
of +1 in all time indices of the state equations but no change 

in the output equation. 
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(3.16b) 



b 3 a n a 22 " a 12 a 21 



c = Y = c 1 b 1 + c 2 b 2 



(3.16c) 



d = <5 



(3.l6d) 




(3.l6e) 



Assuming that the coefficients a, b, c, d and e are given, 
it is necessary to solve equations (3.16) for the entries 
of the matrices of the state and output equations from which 
(3.13) was derived. Any combination of the a.., b. and c n . 
which yields the given values of a, b, c and e is a valid 
realization. But there are eight coefficients in the matrices 
A, B and C which implies that in some realizations there 
could be as many as eight quantization error sources. Further- 
more if the values of a, b, c and e are given to k-bit accur- 
acy, the solutions to equations (3.16) will not, in general, 
yield coefficients to k-bit accuracy. The approach taken 
by Parker and Hess was to force the solutions to be in k-bit 
accuracy by writing (3*l6a, b, c and e) as 




(3.17a) 




(3.17b) 
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(3.17c) 



c = [b 1 c 1 ] Q + 



■Vz’q 



e 



[ a i2°2 b l ]q + [a 21 C l b 2^q " [a il C 2 b 2^q " t a 22 C l b l^q 

(3.l8d) 



where [ • ] indicates quantization of the product within. 
Since there are four equations in eight unknowns, it is 
possible to choose an appropriate set of four unknowns arbi- 
trarily and solve for the other four. Furthermore, the 
solution is simplified by the intuitive realization that 
quantization of product terms can be avoided if, in the 
products of (3.17) > each term except one is unity or one of 
them is zero. Coincidentally, it turns out that this approach 
also results in only four non-zero, non-unity coefficients, 
i.e., the four coefficients remaining after arbitrarily 
choosing four to be unity or zero. This approach thus 
yields the minimum number of quantization error sources as 
well . 

By applying a set of rules for the selection of 

unity and zero coefficients, Parker and Hess obtained eight 

2 

system matrices which would meet the criteria. Only two 
of these were unique, two others being transposes of the 



A ; B 
C i 6 



These system matrices were of the form 
partitioning the matrices yields the 
matrices A, B, C and 6. The phrase "system 
array" will be preferred here for reasons to be shown later. 



i.e.. 
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first two, the other four simply involved the interchange 
of rows or columns in the first four. 

The two unique system arrays derived by Parker and 
Hess were 



S = 
a 



-a 



1 



-b 

0 



c e 



1 - 

0 

d - 



- -(1+a) -(1+a+b) 1 ' 

1 10 

c c+e d - 



(3.18a) 



( 3 . 18 b) 



The S array had been previously discussed by Jackson. The 
transposes of these arrays are then 



r -a 




1 

0 

0 



c 



e 

d 




- -(1+a) 
-(1+a+b) 
1 



1 c - 

1 c+e 
0 d . 



(3.19a) 



(3.19b) 



Jackson showed that the transpose of a system array 
yields a configuration identical to the untransposed config- 
uration except that (1) the direction of signal flow is 
reversed and (2) summing nodes become pickoff nodes and vice 
versa. Jackson further showed that for each unique system 
array and its transpose, twelve configurations can be achieved 
through flow graph manipulations. The conclusion of Parker 
and Hess was then that there were twenty-four canonic 
realizations . 
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The previous analysis was based on equating the 
transfer function of equation ( 3 - 3 ) with the transfer function 
of ( 3 . 13 )* A similar analysis can be carried out by equating 
(3.4) to (3.14). Then 

h(z) = (a- + c ' * f 2 ' 1 ,)z -1 = («■ + 4 c ' Z-1 .) 

1 + az + bz ' 1 + az" + 3z 

= H t (z) (3.20) 

where c', d' , e’, <$ 1 , y' and e' are primed to distinguish 
them from the analysis which precedes, while a, b, a and 3 
are not primed since the characteristic polynomial must be 
the same for all realizations. When the coefficients are 
equated, however, it is seen that the equations are identical 
in form to equations ( 3 . 16 ). Thus following the same set 
of rules for selection of unity and zero coefficients, the 
system arrays which result will have the identical form of 
S and S, of equations (3.18). The only exception is that 

cL D 

d will be replaced by d' , c by c’, and e by e', which means 
that for the same transfer function H(z), the arrays which 
result from ( 3 . 20 ) will differ from those resulting from 
(3.19) numerically. In order to indicate which form of the 
transfer function is used to find the entries in the array, 
a T will be used for those resulting from equation (3.4) 
whereas an S designates the ones found by Parker and Hess, 
thus the rationale for calling the groups S and T comes 
from historical development. The next section will force 
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another change in notation (altering the subscripts, e.g.'a' 
in S & ) to avoid an historical conflict and the resulting 
confusion that might be created. 

In summary, the state equations and associated 
transfer functions are given by 

x(n) 
v(n) 



and 

x(n) = Ax(n-l) + B'u(n-l) 

+ H(z) 

v(n) = C*x(n) + 6'u(n-l) 



H t (z) 



d» + 



+ e ’ z 



e ' z 

* -1 . . - 2 , 

1 + az + bz 



-1 



( 3 . 22 ) 



Ax(n~l) + Bu(n-l) 
Cx(n-l) + 6u(n-l) 



H( z) = H g (z) 



d + 



-1 , -2 
cz + ez 

-1 -2 

1 + az + bz 



-1 



( 3 . 21 ) 



E. CANONIC ARRAYS 

When Parker and Hess wrote their state and output equations, 
it was done in such a way as to make it possible to write 
all three equations as one matrix equation. Thus 

x(n+l) = Ax(n) + Bu(n) (3.23a) 

v(n) = C5c(n) + 6u(n) (3.23b) 



63 



because 



x(n+l) 

v(n) 





a ll 


a i2 


V 


II 


a 21 


a 22 


b 2 




1 

1 

: cT 

l 


C 2 


6 



x(n) 
u (n ) 



(3.24) 



Thus, the arrays S a and of equations ( 3 . 1 8 ) could truly 
be called matrices. Examination of the entries in Table 3-1 
shows that by raising the time indices of the state equations 
of line (e) by 1 to coincide with (3.23a), the same result 
is achieved. It is also possible to write the state and 
output equations of lines (c) or (h) directly as one matrix 
equation. These three lines are all in group S. When one 
tries to do the same thing with the equations in group T 
it is found to be impossible. Therefore, the arrays asso- 
ciated with group T cannot be called matrices. Hence, in 
this work only the description "array" will be used when 
referring to the 3x3 arrays S or T. The transpose of an array 
is defined as though it were a matrix. A, B, C and 6 continue 
to play the role of matrices and will be referred to as such. 

As mentioned earlier, Jackson showed that there were 
six realizations possible for an array and six for its trans- 
pose. This can be accomplished by flow diagram manipulations 
or, equivalently, by manipulations of the equations involved. 
Two of the six are identical except that in one d = 0 and 
in the other, d = 1. In the remaining four it was claimed 



that d had to be unity. In fact. It will be shown that there 
are nine manipulations of the equations possible for the 
state equations of group S and that each allows d to be 
zero or unity. There is a realization for all 36 of these 
cases (the nine manipulations of the array S , nine of S , 
each with d = 0 and d = 1) . Figures 3-1 through 3-18 show 
these realizations. In addition, there are the 36 transpose 
configurations of these. To allow the realization of any 
transfer function with a single realization, composite 
realizations are included which usually require an additional 
summer and which change the configuration from d = 0 to d = 1 
by the use of a single switch. 

Furthermore, the state equations in group T can also 
be manipulated to yield different flow diagrams. There is 
not as much flexibility in this group, however, and only forty- 
four unique configurations have been found (compared to 
72 for group S). Figures 3-19 through 3-22 show the 
untransposed realizations for group T. 

Since previous work has referred to the realizations 
in a manner which is not adaptable to a one-to-one corres- 
pondence with the realizations derived here, the notation 
needs to be changed to avoid confusion. 



switch will be indicated as follows (see Figure 3-23) : 

a) A switch which is normally closed for d = 1 and open 
for d = 0 is represented as a multiplication by d. 

b) A switch which is normally open for d = 1 and closed 
for d = 0 is shown as a multiplication by (1-d). 
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Call an array In the form of equation (3.18a) an 

"M" array j an array in the form of equation (3.18b) is an 

"N" array. When the array is realized for the form of the 

transfer function and state equations of (3.21), precede 

its designation by an "S". In this case the entries of the 

array include the unprimed quantities: c, d and e. Table 

3-2a lists all the SM array manipulations. Table 3-2b lists 

those of the SN array. SMqq compares to the S & array of 

h 

(3.18a). SNgg compares to in ( 3 . 18 b). 

When the array is realized for the form of the trans- 
fer function and state equations of (3.22), precede its 
designation by a "T" . In this case the entries of the array 
include the primed quantities: c', d* and e'. Table 3-3a 
lists all the T'M array manipulations. Table 3-3b shows the 
single array, TN^. This compares to the array with 
primed quantities replacing the unprimed. Similarly TM^ 
compares to with primed entries. 

For group T it will be shown that the transpose 
configurations are not the transposes of Tables (3-3) but 
rather correspond to those of Tables 3-2 by substituting 
primed quantities where applicable. 



Zi 

The subscripting notation will be explained after an 
examination of the process of manipulating the state 
equations . 
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-(b+e) 
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-(1+a) -(1+a+b) 1 1 f -(1+a) -(1+a+b) 1 -(1+a) (c+e)- 
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Table 3-2b Canonic Arrays, SN 



1 

«H O 



1 



T3 



o 



+ 

a$ 



(L> 

+ 

o 



ctf 



l 



o 
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Table 3-3a Canonic Arrays, TM. Table 3-3b Canonic Array, TN 



The development of these realizations follows. (Only 
operations which do not increase the number of non-zero, 
non-unity multiplications are allowed.) 

Figure 3-la shows the flow diagram of the array, 

SM qo with a = 1. The state and output equations for this 
realization are, according to equation (3.21): 

x x (n) = -ax 1 (n-l) - bx 2 (n-l) + u(n-l) (3.25a) 

x 2 (n) = x 1 (n-l) (3.25b) 

v(n) = cx 1 (n-1) + ex 2 (n-l) + du(n-l) (3.25c) 



Figure 3-lb shows the realization of these equations 
for d = 0. This configuration corresponds to the realization, 

I 

S al , of previous work. Figure 3-la corresponds to S &2 . 

Figure 3-lc exhibits the composite realization which 
allows d to be either 0 or 1 . When Parker and Hess said there 
were 24 canonic realizations of a transfer function, the 
truth was, that they found four realizations for a transfer 
function with d = 0 and twenty for a transfer function with 
d = 1. These are never the same transfer function. There 
is a class of functions with d = 0 and a mutually exclusive 
class with d = 1. Correspondingly, there is a class of func- 
tions with d' = 0 and a mutually exclusive class with d' =1. 
Taken together, there can be overlap between the classes of 
d and d* . That is ( if d = 0, d' could be 0 or 1, etc. In 
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total there are four classes given by (d = 0, d' = 0), 

(d = 0, d' = 1) , (d = 1, d» = 0) and (d = 1, d' = 1) which 
include all second order transfer functions. The inclusion 
of a composite realization is in the interest of logistics. 
Each composite realization can be used for any second order 
transfer function. 

Proceeding with the development of the other forms, 
consider equation (3.25c). This equation can be rewritten 
by adding and subtracting the quantity bx 2 (n-l) to yield a 
new output equation 

v(n) = cx^(n-l) + (b+e)x 2 (n-l) + du(n-l) - bx 2 (n-l) 

( 3 . 26 ) 

The quantity -bx 2 (n-l) is one which is formed in calculating 
x^(n) so it is available to be used in (3.26) without adding 
another multiplier to the realization. 

Equation (3.26) together with (3.25a and b) are 
realized in Figure 3-2a for d = 1. In Figure 3-2b, d = 0 
and 3-2c and d show the composite realization. These are 
realizations of SMq^ where the array is 

"-a -b l" 

10 0 

c (b+e)-b d 
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On the other hand it is possible to rewrite equation 
(3.25a) by adding and subtracting ex 2 (n-l) to the right hand 
side. Then 

x.^n-1) = -ax ]L (n-1) - (b+e)x 2 (n-l) + u(n-l) + ex 2 (n-l) 

(3.27) 

Equation 3.27 together with equations (3.25b and c), 
the original second state equation and the original output 
equation, are realized in Figures 3-3a, b, c and d with 
d = 1, d = 0 and two composites, respectively. In this case 
the array is SM Q2 given by 

’-a e-(b+e) 1 

10 0 
c e d 

It is now possible to describe the subscript notation used. 
In SMqq there are no manipulations of the original array. 

In and SMq 2 the second subscript is changed to indicate 

a change in the second column of the array. When the sub- 
script is 1 it indicates adding and subtracting the coeffi- 
cient of the first row in that column to the coefficient 
of the last row in that column. A subscript 2 indicates 
the opposite operation. These operations can be performed 
in the first column as well, or in both columns at the same 
time and in the same or opposite directions. A complete set 
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of arrays with all these operations are given In Tables 3-2a 
and b. Only the untransposed arrays for group S are given 
in the tables. Figures 3-1 through 3-18 show all the untrans- 
posed configurations implied by Table 3-2, whereas Figures 
3-24 through 3-41 show the realizations found from the 
transposes of the arrays of Table 3-2. 

Figure 3-5a is the form historically called the 
"canonic” realization if c*q is factored out of the numerator 
of equation (3.2) . 

For the arrays of group TM, manipulations are only 
possible between the two entries corresponding to a^ and 
c 2 , that is the first entry in column 1 of the array and the 
last entry in column 2. (All manipulations in group S were 
accomplished within the same column.) There are apparently 
no ways to rearrange the equations of array TNq without 
introducing more multipliers and summers. 

To change TMq into TM^ , consider the state and output 
equations for TMq. From equation (3.22), 



x^(n) = -ax^n-l) - bx 2 (n-l) + u(n-l) 


(3.28a) 


x 2 (n) = x 1 (n-l) 


(3.28b) 


v(n) = c'x 1 (n) + e'x 2 (n) + du(n-l) 


(3.28c) 



By adding and subtracting the quantity aXp(n) on the 
right side of (3.23c) the output equation becomes 
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v(n) = c’x 1 (n) + (a+e’)x 2 (n) - ax 2 (n) + du(n-l) (3.29) 

But by equation (3.28b), x 2 (n) = x^n-1) so 
v(n) = c , x 1 (n.) + (a+e')x 2 (n) - ax 1 (n-l) + du(n-l) (3-30) 

The quantity -ax^(n-l) is one which is formed in 
calculating x-^(n) so that it is available to be used in (3-30) 
without adding another multiplier to the realization. 

Equation (3.30) along with (3.28a and b) are realized 
in Figure 3-20 as TM^. The opposite manipulation indicated 
by TM 2 in Table 3-3 is shown in Figure 3-21. 

The untransposed realization of TNq is given in 
Figure 3-22. 

5 

TM^ and TM 2 have no transpose configurations, but 
TMq does. It is not true, however, that the transpose of TMq 
is realized by reversing the flow in Figure 3-19 and inter- 
changing summing nodes and branching nodes. The realization 
is found by applying the state equations of equation (3.22). 
Furthermore, taking Sm!' . , replacing c, e and d by c’, e* and 
d* respectively and using the state equations of (3.22) yields 
another set of realisations for H(z). These realizations are 
designated TM^ • (with i and j corresponding to the subscripts 

1 J 

of SM^\). Note that TMq and TMq Q are the same realization. 

The same correspondence is true for the transpose of 
TNq. Taking SN?j , replacing c, e, and d by c’, e' and d' 

5lf the TMq and TM2 arrays are transposed it is seen that 
five multipliers would be required, thus they would no longer 
be canonic in terms of multipliers and error sources. 
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respectively and using the state equations of 3.22 yields 
the realization TN^j , with TNq Q the same as TNq . 

Figures 3-^2 through 3-59 show all the transpose 
configurations for the T forms. 

F. CONCLUSIONS 

This development has included a very exhaustive search 
of allowable state equations and transfer functions for 
realizations of second order digital filters, canonic in 
the number of delay elements and multipliers (and thus quan- 
tization error sources after multiplication) and also canonic 
in that the multiplier coefficients are linear combinations 
of the coefficients once they are given in the S or T form 
of the transfer function. 

In brief, the arrays of Tables 3-2 and 3-3 when combined 
with the state equations of (3.21) and (3*22), respectively, 
yield the realizations of Figures (3-1) through (3-22). The 
transposes of the arrays of Tables 3-2 when combined with 
both sets of state equations yield the realizations of Figures 
(3-2*0 through (3-**l) for Equations (3.21) and Figures (3-**2) 
through (3-59) for Equations (3.22). 

In total, these realizations consist of 36 realizations 
for transfer functions with d = 1 and 36 for d = 0, as well 
as 22 realizations for d' = 1 and 22 for d' = 0. That is, 
for a particular transfer function, there are 58 realizations 
since d cannot be both 0 and 1 nor can d 1 be both 0 and 1. 
Included are composite realizations which allow d or d' to 
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be either 0 or 1. These require the inclusion of a switch 
to indicate the particular value and sometimes an additional 
summing device. 

The concept of "canonic arrays" introduced here made 
this development possible. If one forces the array to play 
the role of a matrix then only the state equations of the 
S forms are possible. This forces the transfer function 
also to be in the S form. The T forms only arise when one 
allows for the concept of canonic array or, equivalently, 
separating the operations of the state equations from those 
of the output equations. This allows the time index of the 
states on the right hand side of the equations to be different. 
This factor is essential in developing the T forms. 
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Figure 3-la. Realization of SMqq with d-1 




Figure 3-lb. Realization of SMqq with d-0 
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Figure 3-lc. Composite Realization of SMqq with d=0 or 
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• Figure 3-2a. Realization of SM Q1 with d=l 




Figure 3-2b. Realization of SMq^ with d=0 
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Figure 3-2c. Composite Realization of SM fJ j with d=0 or 1 





Figure 3-2d. Composite Realization of SMq-j with d=0 or 1 
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Figure 3- 3a. Realization of SMq 2 with (1=1 




Figure 3-3b. Realization of SMq 2 with d=0 
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Figure 3- 3c. Composite Realization of SM^ with d=0 or 1 




Figure 3-3d. Composite Realization of SM^ with d=0 or 
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Figure 3-4a. Realization of SM^q with d=l 



v(n) 




Figure 3-4b. Realization of SM^q with d=0 
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Figure 3-4c. Composite Realization of SM^q with d=0 or 1 




Figure 3-4d. Composite Realization of SM-jq wtih d=0 or 1 
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Figure 3- 5a. Realization of SM^ with d=l 




Figure 3-5b. Realization of SM-j-j with d=0 
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Figure 3-5c. Composite Realization of SM^ with d=0 or 1 




Figure 3-5d. Composite Realization of SM, ^ wtih d=0 or 1 
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u(n-l ) 




Figure 3-6a. Realization of SM^ with d=l 




Figure 3- 6b . Realization of SM-jo with d=0 
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Figure 3-6d. Composite Realization of wtih d=0 or 1 
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u(n-l) 



*5 






Figure 3-7a. Realization of $M 2n wtih d=l 



u(n-l) 



-(a+c) 







Figure 3-7b. Realization of SM 2 q with d=0 
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Figure 3-7c. Composite Realization of Sf^g with d=0 or 1 




Figure 3-7d. Composite Realization of Sl^g with d=0 or 1 
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u(n-l) 






Figure 3-8a. Realization of Sf-U, with d=l 




Figure 3-8b. Realization of Sf^ with d=0 
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Figure 3- 8c. Composite Realization of Sf^-j with d=0 or 1 




Figure 3-8d. Composite Realization of SM^-j with d=0 or 1 
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Figure 3-9a. Realization of Sf^ with d=l 



u(n-l) v(n) 




Figure 3-9b. Realization of Sf^ with d=0 
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u(n-l) 






e 







-(b+e) 









Figure 3-9c. Composite Realization of with d=0 or 1 




Figure 3-9d. Composite Realization of Sf^ with d=0 or 1 







Figure 3-10a. Realization of SNqq with d=l 




Figure 3-10b. Realization of SNqq with d-0 
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Figure 3-1 Oc. 



Composite realization of SNqq 



with d=0 or 1 
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Figure 3-1 la. Realization of SNq^ with d=l 




Figure 3-1 lb. Realization of SN^ with d=0 
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Figure 3-1 1 c. Composite Realization of SNq-j with d=0 or 1 




Figure 3-1 Id. Composite Realization of SN'qi wtih d-0 or 1 
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u(n-1) v(n) 




Figure 3-1 2a. Realization of SNq,, with d=l 




Figure 3-1 2b. Realization of SN^ with d=0 
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Figure 3-12c. Composite Realization of SNq£ with d=0 or 1 




Figure 3-i2d. Composite Realization of SMg^ with d=0 or 1 



100 





Figure 3-13a. Realization of SN^ with d=l 



v(n) 




Figure 3-1 3b . Realization of SN^q with d=0 
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Figure 3-13c. Composite Realization of SN^ with d=C or 1 
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u(n-l) 







Figure 3-1 4a. Realization of SN^ with d=1 




Figure 3-1 4b . Realization of SN^ with d-0 
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Figure 3-1 4c. Composite Realization of SN^ with d=0 or 1 
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u(n-l) 




Figure 3-1 5a . Realization of SN^ with d=l 




Figure 3- 15b. Realization of SN'^ with d=0 
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Figure 3-1 5c. Composite Realization of SN^ with d=0 or 1 



u(n-l) 






Figure 3-15d. Composite Realization of SN^ with d=0 or 1 
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Figure 3-16a. Realization of S^q with d=l 




Figure 3-1 6b. Realization of SN^q with d=0 
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Figure 3-16c. Composite Realization of S^q with d=0 or 1 




Figure 3-16d. Composite Realization of SN 20 with d=0 or 1 



108 




Figure 3-1 7a. Realization of with d=l 




Figure 3~17b. Realization of SN with d=0 
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Figure 3-17c. Composite Realization of S^, with d=0 or 1 




Figure 3-1 7d. Composite Realization of S^-j with d=0 or 1 
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Figure 3-18a. Realization of SN^ with d=l 




Figure 3-1 8b . Realization of SN^ with c!=0 
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Figure 3- 18c. Composite Realization of SN^ with d=0 or 1 




Figure 3-1 8d. Composite Realization of Sfl 22 with d=0 or 1 
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u(n-l) 




Figure 3-1 9a. Realization of TMq with d ' =1 




Figure 3-1 9b. Realization of TMq with d'-O 
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Figure 3-19c. Composite Realization of TMq wtih d 1 =0 or 1 



11*1 




Figure 3-20a. Realization of TM-j with d ' =1 




Figure 3-20b. Realization of TM-j with d'=0 
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Figure 3-20c. Composite Realization of TM-j with d'=0 or 1 



u(n-l) 






Figure 3-20d. Composite Realization of TM-| with d'-O or 1 




Figure 3-21 a. Realization of Tf^ with d'=l 




Figure 3-21 b. Realization of T with d'=0 
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Figure 3-21 c. Composite Realization of TMg with d'=0 or 1 




Figure 3-21 d. Composite Realization of Tf^ with d 1 =0 or 1 
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Figure 3-22a. Realization of TNq with d'=l 




Figure 3-22b. Realization of TN Q with d'=0 
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c' 



c'+e' 



Composite Realization of TN Q with d'=0 or 1 



Figure 3-22c. 



mul tipi ier 



-G> 



N.C. (d=l ) 



switch 



AND gate 

V „ for d=l 
cc 

6ND for d=0 
(a) POSITIVE LOGIC 




<±[> 



multiplier 



,N.0.(d=l) 



switch 



AND gate 

GND for d=l 

V for d=0 
cc 

(b) NEGATIVE LOGIC 

N.O. = normally open 
N.C. = normal ly closed 

Figure 3-23. Equivalent symbols for digital switches 
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CO 



u(n-l) ( 



v(n) 







Figure 3-24. Realization of SM Q q with d=0 or 1 




Figure 3-25. Realization of SMq^ with d^O or 1 
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Figure 3-26. Realization of SMq^ with d=0 or 1 




123 




u(n-l) 



6 





v(n) 

d=0 > 



Figure 3-28. Realization of with d=0 or 1 




Figure 3-29. Realization of SM^ with d=0 or 1 








Figure 3-31. Realization of Sl^ with d=0 or 1 
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u(n-l) 



Figure 
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Figure 3-33. Realization of SN Q q with d=0 or 1 




Figure 3-34. Realization of SNq^ with d=0 or 1 
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Figure 3-36. Realization of SN^q with d=0 or 1 
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Figure 3-37. Realization of SN^ with d=0 or 1 
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1 




Figure 3-39. Realization of S^q with d=0 or 1 




Figure 3-40. Realization of with d-0 or 1 
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Figure 3-41. Realization of SN^ with d=0 or 1 
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Figure 3-42. Realization of TM Q g with d'=0 or 1 
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Figure 3-46. Realization of TM^ with d'=0 or 1 




Figure 3-47. Realization of TK^ with d‘=0 or 1 
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u(n-l) 

■ — — ■ i 



Figure 




;-50. Realization of Tf^ with d ' =0 or 1 
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Figure 3-51. Realization of TN q q with d ' =0 or 1 
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Figure 3-54. 



Realization of TN^ with d'=0 or 1 
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t 



Figure 3-55. Realization of TN, j with d'=0 or 1 

CD — j. 




v(n) 



Figure 3-56. Realization of TN^ with d'=0 or 1 



139 








Figure 3-57. 



Realization of T^q with c!=0 or 1 
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IV. ERROR ANALYSIS OF RECURSIVE DIGITAL FILTERS 



A. INTRODUCTION 

The effect of finite precision arithmetic on the perfor- 
mance of digital filters has been considered from the deter- 
ministic (limit cycle) and probabilistic (random noise) point 
of view by several authors [3 , 44-48] . Calculation of output 
noise due to quantization, such as round-off after multipli- 
cation of signals by constant factors, usually involves a 
difficult contour integration [48] which does not readily lend 
itself to physical interpretation in terms of filter structure 
or configuration for purposes of comparison. In addition, 
the usual noise formulas are based upon the assumption that 
the quantization error sources within the filter are uncor- 
related. Although this assumption has been confirmed as 
reasonable by several authors [35,50], it is interesting to note 
that under the conditions of a limit cycle, quantization 
errors within the filter are highly correlated. Thus the 
influence of correlation on noise calculations remains a 
question . 

In this chapter, several results are presented. First, 

m 

following a procedure introduced by Parker and Yakowitz [3]> 
but without assuming uncorrelated error sources, the covariance 
matrix of the state errors and the variance of the output 
error are expressed in terms of the correlation matrices of 
the quantization error sources This enables the effects of 

1-This approach was presented at the Arden House Workshop 
on Digital Signal Processing [57]. 
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correlation between sources to be taken into account. The 



results are presented in a general form which enables the 
structure of the filter to be introduced independently. 

More specifically, the expression for the output error 
variance is applied to second order sections and then to the 
canonic realizations of Chapter III. 

Experimental verification of the results are presented 
in Appendix D . 

B. STATE SPACE FORMULATION 

In Chapter III the equations for a digital filter in 
state variable form are given by 



x(n) = 
v(n) = 
for the S forms 
x(n) = 
v(n) = 



Ax(n-l) + Bu(n-l) 

Cx(n-l) + 6u(n-l) 

(see Chapter III, Equations 
Ax(n-l) + B'u(n-l) 

C'x(n) + 6'u(n-l) 



(4.1a) 

(4.1b) 



(3.21)) 2 , and 



(4.2a) 



(4.2b) 



2rp he use of the terminology S forms (or T forms) is not 
restricted to second order filters. The state and output 
equations in matrix form are in general applicable to n 1 
order filters. Even to the transfer function of an n order 
filter, the S,T terminology applies in general. S means 
dividing the numerator by the denominator in ascending powei s 
of j T means to divide in descending order. 



for the T forms (from Equations (3.22)), where A, B, B', C, 

C' , 6 and 6’ are coefficient matrices; x(n) , the vector of 
state variable' at time nT; u(n), the input; and v(n) the 
output . 

If one assumes finite precision arithmetic in carrying 
out the digital filter algorithm, then Equations (4.1) become 



x*(n) = [Ax*(n-1)] + 

- ~ q 


[Bu(n-l)] q 


(4.3a) 


v«(n) = [ Cx* (n-1 ) ] + 


[6u(n-l)] q 


(4.3b) 



where the asterisk denotes rounded signal variables and 

[ ] denotes the process of quantization after multiplication 

Q 

of a previously quantized signal variable by a multiplier. 

% h 

That is, if x.(n) is the i component of x*(n) and a., is 

1 1 J 

the ij th entry in A, then the i th component of [Ax*(n-l)l 

P Q 

is given by E [a..x.(n-l)] where p is the order of the 
j=l 1J J q 

filter. 

Following Yakowitz and Parker [ 2 ], the following bounded 
vectors can be defined: 



e(n) = Ax*(n-1) - [Ax*(n-l)] q + Bu(n-l) - [Bu(n-l)]q 

(4.4a) 



^Algorithms which quantize after summation are also 
possible and have been studied [ ] . Further discussion of 

this and other’ methods are considered In Section D. 



e(n) = Cx*(n-1) - [Cx*(n-1)] + 6u(n-l) - [6u(n-l)] 

- - q q 

(4.4b) 

These are the quantization errors introduced into the filter 
at each time, nT. There can be several contributions to each 
component of e(n) depending upon the number of multipliers 
in the appropriate row of A and B which require the product 
to be quantized; and similarly for e(n) . 

If one defines the error, y(n), as the difference between 
the states of the finite precision realization and those of 
the ideal infinite precision realization, then 

y(n) = x(n) - x* (n) 

= Ax(n-l) + Bu(n-l) - [Ax*(n-l)] q - [Bu(n-l)] q 

(4.5) 

by Equations (4.1a) and (4.3a). The error, y (n) , is by 
definition the accumulated error in the states of the finite 
precision realization. In order to find an equation for the 
propagation of this error. Equation (4.5) can be rewritten as 

* 

y (n) = Ax(n-l) + Bu(n-l) - [Ax*(n-l)] q - [Bu(n-l)] q 
+ Ax*(n-1) - Ax*(n-1) 

= A[x(n-1) - x* (n-1) ] + Ax*(n-1) - [Ax*(n-l)] q 
+ Bu(n-l) - [Bu(n-l) ] 
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y (n) = Ay(n-l) + e(n) 



(4 .6a) 



Thus, the accumulated error in the states at time nT depends 
on the previously accumulated error (at time (n-l)T) and the 
error introduced by quantization at time nT. 

Applying Equation (4.6a) recursively and assuming that 
there is no accumulated error at time nT = 0, i.e. x*(0) = x(0), 
then the dependence of y(n) on all previously introduced 
quantization errors can be found to be 



Equations (4.6a-c) show that the effect of quantization 
errors introduced "long before" time nT becomes small since 
the filter is assumed to be stable and the characteristic 
equation of the filter depends solely on the matrix A. What 
may be considered "long before" depends strongly on the 

A 

characteristic of the filter. In. particular, for second 
order filters, a low damping coefficient would imply a con- 
siderable time for the effect of a quantization error to 
become negligible. This idea will be useful in evaluating 
the results of the derivation and will be seen to be applicable 




(4.6b) 



or by a change of variable in the summation 



y(n) = Z A e(n-i) 
i=0 



(4.6c) 
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to the stochastic lormulation in Section C as well as in 
the deterministic Equation (4.6). 

Defining Av(n) as the difference between the outputs of 
the finite and infinite precision realization, a similar 
derivation yields 

Av(n) = v(n) - v*(n) 

= Cx(n-l) + 6u(n-l) - [Cx*(n-1)] - [5u(n-l)] 

+ Cx*(n-1) - Cx*(n-1) 

= Cy(n-l) + e(n) 
n-1 . 

= C E A n_1 “ J e(j) + e(n) 
j=l 

n-2 . 

= C E A 1 e(n-i-l) + e(n) (4.6d) 

i=0 

The quantity, Av(n) is the accumulated error at the out- 
put of the digital filter at time nT. Equation (4.6d) shows 
that this accumulated error depends on all previous quanti- 
zation errors generated in the state equations (not including 
the errors generated at time nT) but Av(n) only depends on 
the present quantization errors generated in the output 



equation . 



Equations (4.6) apply to. the S forms of the state eaua- 
tlons . In order to do the same for the T forms it is neces- 
sary to redefine the quantization errors in the output 
equation at each time nT as 

e’(n) = C ' x * ( n ) - [c'x*(n)] q + 5' u (n-l) - [6'u(n-l)3 

(4.7) 

Keeping the definition of e(n) the same as (4.11a) except 
to replace B by B' the accumulated state and output error 
equations become 

n n-i n_1 i 

y(n) = Ay(n-i) + e(n) = E A J e(j) = E A e(n-i) 

j=l i=0 

(4.8a) 



and 



Av(n) = C ’y (n) + e ' (n) 



n n , 

= C' E A n " J e(j) + e ’ (n) 

j=l 



= C' E A e(n-i) + e'(n) (4.8b) 

i=0 

Equation (4.8a) is identical to Equations (4.6a-c) but 
(4.8b) has a subtle difference. This difference is that, for 
the T forms, the output error depends on all previously 
generated errors in the state equations including the most 
recent (i.e. at time nT) . The error Av(n) still only depends 
on the most recent error e'(n) of the output equation. 



The difference in the expressions for Av(n) for the S 
and T forms is a result of the choice of state and output 
equations and will present no difficulty in the stochastic 
development which follows. 

C. STOCHASTIC FORMULATION WITH CORRELATION 

A general expression for the variance of the output error 
when the quantization errors are assumed to be discrete ran- 
dom processes can be developed by making use of the error 
propagation equations of the previous section. This expres- 
sion can then be used for analyzing the specific cases of 
the second order filter and particularly the canonic forms 
of Chapter III . 

1 . General Formulation 

Consider first the T forms of the state equations and 
thus Equations (4.8) as expressions for the accumulated 
state and output errors. Assuming for the present that all 
random variables and processes are zero mean the variance of 
the output error can be expressed as 

a A 2 v (n) = E[(Av(n))(Av(n)) t ] (4.9) 

where E [ * ] is the expected value. 

By applying the first equality of Equation (4.8b), 
the output error variance becomes 



°Av (n) = E[(C 'y (n) + e» (n))((e* (n))* + y t (n)(G’) t )] 

= C'A y (n)(C') t + C' E[y(n)(e' (n)) t ] 

+ E[e* (n)y t (n) ](C' + E[e ' (n) ( e ' (n) ) ] 

= C'A y (n)(C') t + 2C ' E[y (n) e ' (n) ] + a* (n) 

(A. 10) 

where 



A y (n) = E[y(n)y t (n)] (4.10a) 

2 2 

a c (n) = E[(e'(n)) ] for the single output case 

(4.10b) 

and 



E[e' (n)y^ (n) ](C' ) t = C 'E[y (n) e ' (n) ] (4.10c) 



Equation (4.10c) arises when a scalar output is 
assumed; thus the cross terms in (4.10) are combined in the 
final expression of (4.10). 

A (n) is the covariance matrix of the accumulated 

y 

state errors at time nT. The dependence of A (n) on the 

kJ 

statistics of the quantization errors is found by substituting 
(4.8a) into (4.10a). Then 
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A (n) = E[{ 2 A n-i e(i)H 2 e fc ( j ) (A t ) n_J } ] 

y 1=1 j=l 

= 22 A n_i E[e(i)e t (j)](A t ) n_J (4.11) 

1=1 j=l 



Defining a quantization error correlation matrix as 



R e (n,m) = E[e(n)e t (m)] (4.12) 



enables (4.11) to be written as 



A (n) =22 A n_i R (i,j)(A t ) n 3 (4.13) 

y i=l j=l 



It is now possible to write the expression 



AA (n)A t =22 A n “ i+1 R (i,j)(A t ) n " J ’ + 1 
y i=l j=l 

n-1 n-1 ^ n i 

= E EAR (i+l,j+l) (A ) 3 (4.14) 

i=0 j = 0 



Subtracting (4.14) from (4.13) yields 



A (n) - AA (n)A t = R (n ,n) - A n R (1,1) (A t ) n (4.15) 

y y e c 

•• n 2 1 R p (n,j)(A t ) n_J ' - n E 1 A n R ( 1 , j+1 ) ( A* ) n “ 
j=l 6 J=1 

+ ^V^R (i,n) - n E 1 A n " i R (i+l,l)(A t ) n 
i=l e i=l 

+ V n 2 1 A n " i [R p (lJ) - R (i+l,j+l)](A t ) n 
i=l j=l 6 
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Equation (4.15) can be simplified if it is assumed 



that all the processes involving quantization errors are wide 
sense stationary and mutually wide sense stationary. Under 
this assumption, the correlation matrix R e (n,m) is a function 
only of the difference of its arguments, (m-n) (taken for 
notation purposes in that order) . Then the correlation 
matrix can be written 



the first and second terms of (4.15) contain R e (0). Further- 
more, the double summation term becomes the 0 matrix, so that 
(4.15) reduces to 



R (n,m) = R (m-n) 



(4.16) 



Under the assumption of wide sense stationarity 



A (n) - AA y (n)A fc = R 0 (O) - A n R 0 (O)(A t ) n 



+ Z R ( j -n) ( A 



n-l/v 




n-1 

+ Z A 
i=l 




n " i R e (-i)(A t ) n 

(4.17) 



Now consider the entries of the correlation matrix. 

If e . (n) and e.(m) are the i th and j th components of e(n) 
i J 

and e (m) respectively then the ij entry of the correlation 



matrix is 
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(4.18) 



[R e (n J m)] i j = [E[e(n)e t (m) ] ] 

= e i (n)e^(m) = e^CmJe^^Cn) 

= [E[e(m)e t (n)]]., = (X(m,n)]., 

J -L “ J 1 

This implies that the transpose of the correlation matrix 
is the correlation matrix with its arguments reversed. For 
the case of mutually wide sense stationary processes this 
means 



H e OO = Re ( " k) 



(4.19) 



Making use of (4.19) In (4.17) it can be written 



that 



A (n) - AA (n)A t 

y y 



= R e (0) - A n R e (0)(A t ) n (4.20) 

+ n z 1 {A n " i R (i (n-i) + CA n " 1 R o (n-i)] t } 
i=l e 



n “l v, ^ +. « 

- I {A n ” 1 R (-i) (A; 11 
i=l 



+ [A n " i R e (-i)(A t ) n ] t } 



Now as n becomes large the last summation of (4.20) 

becomes negligible since the filter is assumed to be stable 

which implies that lim A n = The second term in (4.20) 

n-*-°° 

also becomes negligible so that in the limit 
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(4.21a) 




AA A fc = R (0) + L { A 1 R ( 1 ) + [ A^R (i)]*} 
y i=l e e 



or 



A y - AA y A^ = R e (0) + ^^{A 1 R e (D + R e (-i)(A t ) 1 } (4.21b) 

The left hand side of Equations (4.21) can be solved 
for A so that these equations essentially give an expression 
for A , the covariance matrix of state errors, in terms of 

y 

the correlation matrices of the quantization error sources 
for large n under the assumptions of zero-mean, wide sense 
stationary and mutually wide sense stationary error processes. 
The solution to 

A - AA A t = F (4.22) 

y y 

(where F is a symmetric matrix which represents the right 
hand side of (4.21)) can then be substituted into (4.10) as 
part of the solution for the output error variance. 

The factor E[y(n)e'(n)] which also appears in (4.10) 
must be found before the output error variance can be written 
solely in terms of the statistics of the quantization error 
sources. From (4.8a) it follows that 

n-l , 

E[y (n)e ' (n) ] = I A 1 E[e(n-i)e ' (n) ] (4.23) 

i=0 
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Now define a cross correlation matrix (in this case 



a vector) between £(n) and e'(m) as: 
R ,(n,m) - E[e (n) e ' (m) ] 



.2it ) 



Under the assumption of wide sense stationarity , as 
before, the cross correlation vector depends only on the 
difference in its arguments, i.e. 

R , (n,m) = R , ( m-n ) (4.25) 

ee ' ee 



so that (4.23) becomes 



E[y (n)e’ (n) ] 



n-1 , ^ 
l An i (i) 
i=0 ee 



and it follows that 



(4.26) 



lim E[y (n) e' (n) ] 
n-n» 



Z A R 
i=0 



ee 



(i) 



(4.27) 



Substituting (4.27) and the solution to (4.22) into 
(4.10) yields the output error variance for the T forms. 



°Av ( £> e,) = 



C 1 A (e) (O') 

O' 



+ 20 ’ Z 
1=0 



A 1 ^ i ( i) 

ee 



+ o 



2 (4.28) 



as a function of the statistics of e_ and e' . 

For the S forms, the difference in the output error 
propagation equation (4.7c) causes a slight difference in the 
expression for the output error variance . 
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Substituting (4.7c) into (4.9) yields 



o A 2 v (n) = E[(Cy(n-l) + e(n))(e t (n) + y^n-Dc*)] 

= CA (n-l)C t + 2C E[y (n-l)e (n) ] + o 2 (4.29) 

y ^ 

The first term of (4.29) contains A (n-1) instead of A (n) 

y y 

which appeared in (4.10). But in the limit 

A (n-1) = A (n) = A (e) . The second term of (4.29), however, 

y y y 

results in a slightly different expression in the summation 
form. Using (4.7b) 

n-2 

E[y (n-l)e (n) ] = Z A E[e(n-l-i)e(n) ] 

i=0 

n-2 j ~ 

= Z A R (i+1) (4.30) 

i=0 ee 

which in the limit becomes 

°o 

lim E[y (n-1) e(n) ] = ZA i R (a (i+l) (4.31) 

n+~ ~ i=0 e 

Then the expression for the output error variance for the 
S forms is 

o a 2 v (e,e) - CA^eJC 1 + 2C Z^A^^i+l) + 0 2 0.32) 

Equations (4.28) and (4.32) together with (4.21) are 
the principal results of this section. They express the output 
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error variance of the S and T forms in terms of the statistics 
of the internally generated quantization errors including the 
effects of correlation. The assumptions made were that all 
processes had zero mean and were wide sense stationary and 
mutually wide sense stationary. 



the expressions do not change appreciably. All that is 
necessary is to redefine the covariance and correlation 
matrices : 



When .the assumption of zero mean processes is relaxed 



A (n) = E[{y(n) - m(y)Hy(n) - m(y)} t ] 

J 



E[y(n)^_ t (n)] - m(y (n) )v^ (y (n) ) 



A (n) - m(y(n) )m t (y (n) ) 

o 



(^ . 33 ) 



R (n,m) - E[{e(n) - m(e(n) ) } {e (m) - mCeCm))} 1 ^ 



= E[e(n)e t (m)] - m(e(n) )m t (e(m) ) 




( 4 . 34 ) 



R 



(n,m) - E[{e(n) - m(e (n ) ) } {e (m) - y(e( m))}] 



= E[e(n)e(m)] - m(e (n) ) y ( c (m) ) 



R (n.m) - m(e(n) )y(e(m) ) 
ee — — 



( 4 . 35 ) 
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where m and y are the means of their arguments. 

When it is assumed that the processes are wide sense 
stationary it is possible to write 



where the dependence on a time index has been removed because 
wide sense stationarity implies dependence- only on time 
differences . 

Then the expressions for the output error variance 
for the S forms is 



Ay (e) = A y (e) - m(y)m t (y) 



(4.3 6) 



R e (n,m) = R e (m-n) = R e (m-n) - m(e )m t (e_) 



(4.37) 



R ee (n,m) = R ee (m-n) = R 0e (m-n) - m(e)y(e) (4.38) 




(4.39) 



and for the T forms 



• A 

•l- 2C' Z A 1 R ee i (i) + 0 



00 



a A v ( - ,e,) = C'Ay(e)(C') t 



(4.40) 



where A (e) is the solution to 
«y 



A (e ) - AA r (e)A t = R p (0) + Z {A 1 ^ (i) + [A 1 ^ (D^J 

y - y - e 3-i u 



(4.41) 
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and 



o 2 z = E[ ( £ - y £ ) 2 ] = E[e 2 ] - y 2 (4.42) 

These equations ((4.39) through (4.42)) give the 
general expression for the output error variance assuming 
only wide sense stationarity . The effects of structure are 

£ A Q 

imbedded in the nature of R (i), R (i) and c and in the 
entries of the matrices of the state and output equations. 

The next section examines these effects. 

a. Structural Effects (General) 

In order to take topological or structural effects 
into account, it is necessary to consider the terms of the 

correlation matrices cf error sources as given by (4.16) and 

4 * 

(4.25). When zero correlation is assumed, R g (k) is taken 

A 

to be zero for all k / 0, and R (0) is taken to be a diagonal 

c 

matrix where the k entry on the diagonal is the sum of the 
second moments (variance) of each individual error source 
comprising e^(n) . Assuming uniform distributions, these 

c 

terms become-' (o ) times the number of unique non-zero non- 
unity multipliers in the k 1 row of the matrices which 
contribute to the error e^(n) . 

When correlation is assumed (both spatial and 
time) the picture becomes much more complicated. In general, 

^This development returns to the assumption of zero mean 
processes for ease of notation since there is no loss of 
generality in doing so. 

5 c 2 = h 2 /12 for rounding. h 2 /3 for truncation. 
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the (ij) th 



term of the correlation matrix is 



given by 



R (n,m) - E[e(n)e t (m) ]1 _ r . . 

L e Jij L “ - J ±j = ELe 1 (n)e ( 



m)] 



= E 



Pi Qj 



k 



1 Z e. (n)e.n(m) 
=1 £=i ^ J* 



P i 

- Z L E[e lk (n)e J . p (m) ] 



k=i a=i 



(4.43) 



where p ± is the number of multiplier errors contributing to 
Qj is the number of multiplier errors contributing to 
e.; and k,£ are additional indices to identify the contribu- 
tions to e. , e. respectively; that is 

-L J 

P i Q j 

e.(n) = Z e., (n) and e,(n) = Z e. f (n) (4.44) 
1 k=l 1K J l=i ^ 



In all, then, there are six indices to be considered to 
determine the entries of the covariance matrix; two (i,j) 
associated with states; two (k,&) associated with multipliers 
and two (n,m) associated with time. When wide sense station- 
arity is assumed, then one index representing time difference 
replaces the two indices of time. 

Considering a typical contribution to R (n,m), 
such as E[e ik (n)ej ^ (m) ] , the question to be asked is when is 
this correlation statistic non-zero? One obvious answer is, 
when the multiplier coefficients are the same number and the 
signal quantities are identical (e.g., x^(n) = x^Cn-l)). In 
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this case there is a unity correlation coefficient and 

E[e ik (n)e j* (m)] = E[e. 2 k (n)]E[e^(m)] . (4. 45) 

Another possibility to be investigated is when 
the signal quantities are identical, but the multiplier 
coefficients are different. If the signal quantities have 
Gaussian distributions then the correlation coefficient is 
a function of the ratio of the two coefficients as shown in 
Appendix B . As expected in this case if the multiplier 
coefficients are the same, the error statistics have unity 
correlation. Note that the correlation coefficient depends 
on the sign of the ratio of the multipliers. 

Y/hen the multiplier coefficients are the same 
but the signal quantities are not identical, then it is 
suggested that the correlation of errors will depend mainly 
on the correlation between the signal quantities. 

Finally, the most complicated case is when the 
coefficients are different and the signal quantities are 
not identical. It is expected in this case that the 
correlation is small. 

The foregoing discussion applies equally well 
to the terms comprising R (n,m) . Experimental determina- 
tion of correlation effects are considered in Appendix D 
for second order canonic filters . In the next section the 
general error variance equations are simplified for the 
second order filter. 
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2. Simplification tc the General Second Order Section 



a. For Arbitrary State Equations 

For the second order filter the solution to (4.22) 
is a symmetric matrix 

X 

^12 

(4.46) 



A = 

y 



y ll y 12 



22 



where 



Xy i:L - {f 11 [(l-a 11 a 22 ) (1 a 22 ) - a 12 a 21^ 1+a 22 ^ 



+ 2f 12 ta i2 a 2i a 22 + a li a !2 (1 a 22 ) ‘ ] 



(4.46a) 



+ f 22 t ( 1+a 11 a 22 a i2 a 21^ a i2^ /A 



'12 



Xy 2i {f il^ 1-a 22 )a H a 21 + a l2 a 21 a 22 ^ 



+ f 12 C (l-a 1 ^) (l-a 22 ) - a 12 a 2 j] 



+ f 22 [(l“ a 11 ) a 12 a 22 + a ll a 12 a 21' ] J/A 



(4.46b) 



Xy 22 = {f ll [(1+a H a 22 a 12 a 21 )a 21 ] 



+ 2f 12 t a 11 a 1 2 a 21 + a 21 a 22 ^ 1 a ll^ 



(4 .46c) 



+ f 22 [ ( 1_a 11 a 2 2 ^ ^ 1_a ll^ a 1 2 a 21 (l +a 11 ) 
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/ 



with 



A (l “ a n ) ( 1 “ a 22 ) (1-a n a 22 “ a !2 a 21 ) 



2(1-a il )a i2 a 2i a 22 



2(1 a 22^ a H a i2 a 21 



^ 1-a H a 22“ a l2 a 21 )a i2 a 21 



Z|a li a i2 a 2i a 22 



(4.46d) 



and 



F 



f 

f 



f 

11 12 

f 

21 22 



with f 




(4 .H6e) 



For general second order state equations this 
solution could then be substituted into (4.28) or (4.32). 
b. For Canonic Arrays 

When the state equations are specifically those 
of the canonic arrays in Chapter III, it is possible to 
substitute the coefficients found there into Equations (4.46) 
to find A in terms of those coefficients and the entries of 

y 

the matrix F. If the factor E[y_(n-1) e (n) ] in (4.29) is 
given by 

00 

* E[y (n-l)e(n) ] = l A^Ji+l) (4.47) 

1 = 0 
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then the expression for the output error variance for the S 



forms of the state equations become: 



SM: °Av = ^ f u^ c2 + e 2 )(l + b) + ce ( -2a) ] 

+ 2 f 12 [(c 2 + e 2 ) ( ab ) + ce(l-a 2 -b 2 )] 

+ f 22 [c 2 b 2 (l+b) - 2ce(ab 2 ) + e 2 [(l+b) - a 2 ( 

+ 2(cg 1 + eg 2 ) + a * (4 

SMS a A 2 y = {f 1;L (l+b) - 2af 12 + f 22 (l+b)}/A (4 

SN : o ^ ^ — {f^-^CCc t e )(l+b) + cc(—2a)3 

+ 2f 12 [(c 2 + e 2 ) (1+b+ab) + ce ( l-a 2 -b 2 -2a) ] 

+ f 22 [c 2 {2(l+b+ab)-(i+b) 2 (l-b) } 

+ 2ce{(l-b) ( 1+b+ab ) -a (1+a+b) } 

+ e 2 {2(l + b + ab) - a 2 (l-b)}]}/A 

+ 2{c(g 1 + g 2 ) + e(g 2 )} + a 2 (4 

SN t : a A 2 y = {2f 11 d+a+b) - 2f 12 (l+a+b) + f 22 (l+b)}/A 

('I 



1-b ) ] } /A 

.48a) 

.48b) 



.48c) 

.48b) 
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where 



A = (l-b)[(l+b) 2 - a 2 ] (4.48e) 

When the T forms of the state equations are used, 
let the factor E[y(n)e'(n)] in (4.10) be given by 

00 

= E[y (n)e ’ (n) ] = Z A 1 ^ .(i) (4.49) 

i=0 ee 

For the T forms, the expressions for the output error variance 

are identical to Equations (4.48) for the coresponding canonic 

arrays (e.g., TM corresponds to SM , TM to SM etc.) 

/ 2 

except that c, e, g^ g 2 and c £ are replaced by c’, e‘, g^, 
g 2 ' ana a £ 2 , , respectively. 

Note that the transpose forms do not apparently 
depend on c and e. The factors f^. however do depend on 
these coefficients since in the transpose case they appear 
in the state equations rather than the output equation. 

The quantization errors they cause will appear in e(n) thus 

A 

in R (i) and thus in f... 

G 1J 

The output error variance for the transpose cases 

2 

does not however depend on any G or 0 £ since e(n) (ore’ (n)) 
is zero for the transpose case, there being no non-zero, 
non-unity multipliers in the output equation. 

The expressions for the output error variance 
given in Equations (4.48) are identified without subscripts 
but rather by an identification which implies the applicability 



G * = 



g l* 

So' 
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of that expression to all members of the set of realizations 



with that basic identification (e.g. SN 1Q is a member of 
the set of SN realizations). The reason these expressions 
apply to the entire set is that the development to this point 
has depended solely on the basic state equations. The dif- 
ferences of structure implied by the difference between the 
eight basic forms of the state and output equations exhibit 
themselves in the four basic expressions of 0 08) for the 
S forms and the same four expressions for the T forms. The 
differences of structure imposed by the various manipulations 
of the flow diagrams are entirely embodied in the factors 
f . . and g. . What remains to be examined is how that structure 
exhibits itself in those factors. Structural effects were 
discussed in general in Section C2. Now they will be examined 
in the context of the canonic arrays starting with the simpli- 
fying assumption of uncorrelated error sources. This condition 
will then be relaxed to achieve some better insight into the 
effects of correlation. 

c. Uncorrelated Case 

When all the error sources are assumed to be 
uncorrelated, the right hand side of Equation 0.22) becomes 



F = 



E[e2(n) ] 



0 



0.50) 



0 



ECe^Cn) ] 
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2 

where E[e^(n)] is the sum of the variances of the individual 
multiplier errors contributing to e^n) . 

For example, in the realization of the array, 
SM 2 q, the first state equation is 

x 1 (n) = cCx^n-l) + u(n-l)) - (a+c)x 1 (n-l) + x 2 (n-l) 

(4.51) 

therefore 



e^(n) = c(x*(n-l) + u(n-l)) - [c(x*(n-l) + uCn-l)}]^ 
+ {-(a+c) x*(n-l) ) - [(-(a+c)x*(n-l) )] 0 
+ x|(n-l) - [x|(n-l)] q 

= e c (n) + e (a+c) (n) + 0 (4.52) 



and 



2 

E[e*(n)] = E[(e c (n) + e (a+c) (n)} ] 

= E[e^(n) ] + E C e ( a+c )( n) ^ (4.53) 

since E[e (n)e, . *(n)] is assumed to be zero because the 

c (a+c; 

error sources are assumed to be uncorrelated. 



167 



If one assumes that all error sources are 

a p 

identically distributed with variance a , then for the 
example SM 2Q the variance of e^Cn) is 



E[e*(n)] = 2 o 2 



( 4 . 54 ) 



In general, for the uncorrelated case 



F = 




0 



Ap 

L ° y 2° _ 

where (i = 1,2) is the number 

r.. 

non-unity multipliers in the i n 
For the example, the array is 

c-(a+c) 1 

-b 0 

1 0 




(4.55) 



of unique non-zero or 
row of the canonic array. 



c 

e 




(4.56) 



It appears at first that there are three non-zero, 
non-unity multipliers in the first row. But by referring 
to the flow diagram of this realization (Figure 3-20) or 
Equation (4.51), it is seen that there is only one multiplier 
with c as the coefficient and one with a+c. Thus there are 
two quantizations required, one after each of these two 
multiplications . 
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SM 



In the array 

c-(a+c) -b 1 

1 0 0 

c e d 



20 



(4.57) 



there are three unique non-zero, non-unity coefficients in 
the first row and none in the second. Then for SM 20 

:2 

u 

(4.58) 



F = 



3o' 

0 



0 

0 



whereas for SM 



20 



2o‘ 



0 



F = 



0 2o‘ 



Therefore in counting coefficients for forming F, only 
unique non-zero, non-unity coefficients in each. row are 
included, because the same reasoning is true in all the 
realizations . 

To find the entries of the vector G or G 1 , it 
is necessary to examine the relationships between the state 
equations and the output equation. 

From Equation (4.47) and the definition of R e£ (k), 

oo 

G = l A 1 E[e(n-i-l)e(n) ] (4.60) 

i=0 
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and from (4.49) 



G' = 2 A i E[e^(n-i)e' (n) ] (4 . 61 ) 

1=0 

At first glance. Equations (4 .60) and (4.6l) 
would appear to be zero when error sources are assumed to 
be uncorrelated. But it is to be remembered that e(n) and 
e(n) are not error sources. Rather they are composite errors 
made up of the errors from possibly several sources. In 
fact, in some of the realizations of Chapter III a particular 
source of error may find its way into a component of e(n) 
and also into e(n). 

Consider for example the flow diagram for the 
realisation TM^(see Figures 3-20 ) . The first state equation 

for this realization is 

x 1 (n) = -ax 1 (n-l) - bx 2 (n-l) + u(n-l) (4.62a) 

while the output equation is 

v(n) = c'x 1 (n) - ax 2 (n) + (a+e')x 2 (n) + d'u(n-l) 

= c'x^(n) - axj(n-l) + (a+e')x 2 (n) + d'u(n-l) 

(4.62b) 

The error e- L (n) in this case represents the sum 
of the two errors from quantizing ax^(n-l) and bx 2 (n~l) while 
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e’Cn) is made up of the errors from quantizing c'x^n), 
(a+e’)x 2 (n) and ax^ (n-1) , so that 

e l (n) = e a (n) + e b (n) (4.63a) 



and 



e '(n) = e c ,(n) + e (a+e , } (n) + e a (n) (4.63b) 



Then 



E[e 1 (n)e ' (n) ] = E[e a (n)e c , (n) ] + E[e a (n)e (a+e , ^ (n) ] 

+ S[e a (n)] + E[e b (n)e c , (n) ] 

+ E[e b (n)e^ a+e , ^ (n)] + E[e b (n)e a (n) ] 

(4.64a) 

Under the assumption of uncorrelated error sources all these 

2 

terms are zero except E[e.(n)] so that 

cl 

E[e 1 (n)e'(n)] = E[e|(n) D ? 0 (4.64b) 

This implies that in Equation (4 .61) the first term of the 
summation is not zero. The other terms are zero, however, 
in the uncorrelated case. 

Upon examination of all realizations it turns 
out that only TM^ and TM^ fit this pattern. All transpose 
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configurations require that e(n)(or e’(n)) be zero since 
there are no non-zero, non-unity multiplications in the output 
equation. For the SM and SN forms, although in most cases 
there are coefficients which contribute to both e(n) and 
e(n), the term E[e(n)e(n)] does not appear in Equation (4.60). 
TM q and TN q have no coefficients which contribute to both 
e(n) and e ' (n) . 

Again assuming identically distributed error 

*2 

sources with variance o , then for the realizations TM 1 and 

tm 2 



1 




r : 2 i 


g l 




0 


* 






. S 2 _ 




_ 0 _ 



Although there are only two cases when G’ is 
non-zero for uncorrelated error sources, the exercise in 
finding them will be helpful in understanding the effect of 
structure when correlation is considered. 

Finally, the last term of (4.28) or (4.32) 
involves the variance of the quantization errors contributed 
by the output equation alone, i.e., (or a £ ,). When 
uncorrelated error sources are assumed, this term is the 
sum of the variances of the individual errors. Assuming 
identically distributed error sources. 



a 



2 

e 



= vo 



(4.65a) 
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where v is the number of non-zero, non-unity elements in 
the output row of the array. 

Using the previous results. Equations (4.48) for 
the S forms become: 

° ^ (Vi^/A) [ (c 2 +e 2 ) (1+b) - 2ace] + v}a 2 

(4.66a) 

°Av = (y i + ' J 2 )(1+b)cj2/A (4.66b) 

o a 2 v = (2y 1 (l+a+b) + y 2 (l+b)}a 2 /A (4.66c) 

where A is defined in (4.48e). 

For the T forms these equations become: 

a Av = { (y 1 /A)[(( c ' ) 2 + (e ’ ) 2 ) (l+b)-2ac , e t ]+v}a 2 

(4.66d) 

a Av = ' f (> J 1 /A)[((c') 2 +(e , ) 2 )(l+b)-2ac , e , ]+2c , +v}c 2 

(4 . 66e) 

a Av = (y l + V (1+b)a2/A (4.66f) 

° kv ~ ^2y^(l+a+b) + (1+b ) )a 2 /A (4.66g) 

The fact that y 2 = 0 for non-transpose configura- 
tions is already included in Equations (4.66). 

Table 4-1 lists the results of the analysis 
above for each canonic array for the uncorrelated, identically 
distributed case. 



TMqjTNq : 

TM^ , TM 2 : 
TM b : 
TN b : 



SM,SN: 

SM b : 
SN b : 
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ARRAY 


M 1 


M 2 


V 




a jf 
Av 


formula 

§ 


SM 00’ SN 00 


2 


0 


2 




( 2p+2)a 2 


1 


SM 01 ,SM 10, SN 01 ,SN 10 


2 


0 


3 




( 2p+3)a 2 


2 


SM 02> SM 20’ SN 02’ SN 20 


3 


0 


2 




(3P+2)a 2 


3 


SM 12 y ^^21 3 ^*^1 2 3 ^^21 


3 


0 


3 




(3P+3)° 2 


4 


sM ii ,sN n 


2 


0 


4 




( 2p+4 ) a 2 


5 


SM 22 jSN 22 


4 


0 


2 




(4p+2)a 2 


6 


SM iV™w 


2 


2 


0 




4 (1+b )a 2 /A 


7 


SN 1 Wj ' 


2 


2 


0 




[6(l+b)+4a]a 2 /A 


8 


tm 0 ,tn o 


2 


0 


2 




( 2p ' +2) a 2 


9 


TM X 


2 


0 


3 




(2p’+2c'+3)a 2 


10 


TM 2 


3 


0 


2 




(3p'+2c'+2)a 2 


11 


p = [(c 2 • 


. e 2 


)(1 


+ b) 


- 2ace]/A 




p' = c(cc: 


) 2 + 


(e ' 


' ) 2 )(1 + b) - 2ac'e']/A 




A = (1 


- b)[(l + 


b) 2 


- a 2 ] 





TABLE 4-1 FORMULAE FOR a 2 FOR UNCORRELATED, IDENTICALLY 

DISTRIBUTED, ERROR SOURCES (MEAN-0 ;VARIANCE=a ) 
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Appendix C verifies that none of the formulae 
uncorrelated case violate the requirement that they 
a non-negative number for o* Briefly, the results 
appendix are that for a stable filter 



A > 0 




(4.67a) 


P (or p * ) 


> 0 


(4.67b) 


p ’ + c ' > 


1 

" IT 


(4.67c) 


(1+b) /A > 


l 


(4.67d) 



Examination of the relative magnitudes of the 
? 

expressions for o^ v shows that the smallest value for forrnu- 

^ 2 

lae 1 through 6 is (2p + 2)a and from ( ^ . 6 7b ) the smallest 
this can be is 2a . Formula 7 is at least 4 a * by (4.67d) 

A p 

and formula 8 is at least 2a since a is at least -(1+b) . 

^2 

Formulae 9 through 11 predict a minimum variance of 2a , 

a 2 

2.5 a and 1.5o^ respectively. 

It has been stated empirically [ ] that, of 

the canonic realizations derived by Parker and Hess, noise 
generation is less in the S, (SN) forms than in S (SM), but 

U cl 

that the transpose forms had still less, S^(SN ) being 

+- 4- ^ 

better than S (SM ), i.e. the least error occurs in SN . 
Table 4-1 does not predict this, nor does it contradict it. 
But formula 1 predicts that SMqq and SNqq yield the lowest 



for the 
predict 
of that 
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noise for the untransposed S forms, making no distinction 
between SM and SN. Formula 8 says any of the SN t configurations 
are probably as good as SM Q0 and SN 0Q . Furthermore, that 
same formula indicates that the TN^ forms are comparable 
and formula 9 adds TMq and TNq to the group with a minimum 
noise of 2 a . 

The most intriguing suggestion of Table 4-1 is 
that TM ^ may yield the lowest noise , having a minimum variance 
of 1.5 a . 

But all this is highly conjectural especially 
since only minimum values are being compared. There are 
few sets of coefficients which will result in these minima 
being achieved. For example, (l+b)/A is only minimum at 
the origin in the ( a, b) -plane. The factors p and p ' can be 
minimum for any value of a .and b as long as c and e (or c' 
and e') are both zero which is a degenerate case. 

Furthermore, this analysis was based on the 
assumption of identically distributed, uncorrelated error 
sources. More useful and interesting results can be achieved 
by relaxing these conditions. Furthermore it is closer to 
reality as will be seen in the next section, 
d. Correlated Case 

When the quantization error sources are assumed 
to be correlated, the quantities f . . in Equations (4.48) 
are given by the matrix equation 
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CO 



(4.68) 



F = 



R e (0) 



i=l 



{A 1 R e (i) 



[A i R e (i)] t > 



The entries of R (i) were discussed in general 
in Section C-l-a of this chapter. It was shown there that 
these entries consisted of auto- and crosscorrelation terms 
of the error sources. 

For example, consider the realization TM^q for 
which the array is 

-a 1 c * 

-b 0 e » 

1 0 d' 

and the structure is that of Figure 3-42. The error vectors 
for this realization are given by. 



1 

CD 

M 

'S' 

! 




e a (n) + e c t 


e 2 (n) 




e b + e e t 



(4.69a) 



e(n) = 0 



(4.69b) 



where e (n) is the error introduced by the multiplication 
a 

of x^(n-l) by a and the subsequent quantization of that pro- 
duct at time nT. The errors e^(n), e c ,(n) and e e , (n) are 
similarly defined. Figure 4-la shows a schematic represen- 
tation of the error sources for TMqq. Figure 4-lb shows the 
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composite errors e^(n) and e^in) as composite error sources 
where e^n) and e 2 (n) are given by (4.69a). 

A 

The correlation matrix R e (i) is given by 



R e (i) = E[e_(n-i )e_ fc (n) ] 



E[e 1 (n-i)e 1 (n) ] ECe^n-De^n) ] 
E[e 2 (n-i)e 1 (n)] E[e 2 (n-i )e 2 (n) ] 

(4.70) 



Substituting (4.69) into (4.70) yields 



R e (D 



ll(l) 


R !2 (1) 


21 (l) 


R 22 (l) 



(4.71) 



with 



Rj k (i) - E[ej (n-i)e k (n) ] 

R ll (i) = R aa (i) + R ac’ (i) + R c’a (1) + R c’c’ (i) 
R 12 (i) = R ab (i) + R ae’ (i) + R c'b (i) + R c’e’ (i) 
R 21 (i) = R ba (1) + R bc' (i) + R e’a (i) + R e'c’ (i) 
R 22 (1) = R bb (i) + R be' (i) + R e'b (i) + R e'e’ (i) 



(4.71a) 
(4.71b) 
(4.71c) 
(4.71d) 
( 4 . 71e ) 



where 



(4.71f) 



R aB (jL) " E t e a ( n - 1 )e B (n)] 

under the assumption of wide sense stationarity . 

Appendix D describes an experiment which was 
conducted to check the validity of Equations (4.48) for the 
realization TMqq. The "actual" output error variance for a 
white noise input was calculated by summing the square of 
the difference between the outputs of a double-precision 
realization and one which rounded products to one decimal 
point. (This, of course, required the assumption that the 
process was ergodic (or stationary in a strict) sense in order 
for the time average to be equal to the ensemble average.) 
During the process of rounding products, the error from each 
quantization was saved, correlations were calculated and 

A 

substituted into R g (i) and F was computed. The "predicted" 
output error variance was found using (4.48b) and the a 
posteriori information about the correlation of the error 
sources. Excellent correspondence between the actual and 
predicted values was observed for three sets of transfer 
function coefficients. 

Of course F w as not computed as an infinite sum, 
but a sufficient number of terms were used to observe the 
convergence. Of particular interest is the observation that 
a higher damping coefficient caused more rapid convergence. 
This is to be expected since it is the nature of a stable- 
filter that as the power of the matrix A becomes larger the 
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entries converge to zero and that the rate of this conver- 
gence depends on the damping coefficient. But each term of 
F also includes a correlation matrix. The experiment shows 
that convergence to within 0.1 % is faster when the filter 
is driven by an uncorrelated input than when driven by a 
highly correlated one such as a step (Heaviside) function. 
This too is to be expected since a highly correlated input 

implies higher correlation in the states and thus higher 

\ 

correlation in the error sources. This tends to make the 
entries of the correlation matrices larger, hence a slower 
convergence of F is expected. An heuristic examination 
of this idea follows. 

Consider the nature of the quantization error 



source. Figure 4—2a shows the function Cf(n) J 



->o /m ivi ^ n v> r* 



where f(n) is the number to be rounded and h is the quanti- 
zation step size. The dotted line represents the identity 
function f(n) = f(n). Figure 4.2b shows the error 

e[f(n)] = f(n) - [f(n)] . Note that the error can be written 

Q 

as 



e(n) = e ( f ( r. ) ) = f(n) - hk(n) 



(4.72) 



with 



k(n) = C(f(n) + |)/h] x 



(4.73) 



where [.] means the largest integer less than (•). Equation 
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(4.73) implies that 



(f(n) - |)/h < k(n) < (f(n) + |)/h (4.74) 

which implies that 

hk(n) - ^ < f (n) <_ hk(n) + ^ (4.75) 

Now consider the autocorrelation of the error 

E[e(n-i)e(n) ] = E[(f (n-i)-hk(n-i) ) (f (n)-hk(n) ) ] 

= E[f (n-i)f (n) ] - hE[k(n-i)f (n) + f(n-i)k(n)] 
+ h^E[k(n-i)k(n) ] 

= E[f (n-i) f (n) ] - hE[f (n-i)k(n) ] 

- h(E[k(n-i)f (n) ] - hE[k(n-i)k(n) ]} (4.76) 

If the discrete probability density function 
(p.d.f.) of the random variable f(n) is given by Pp(f(n)) = 

Pr(F = f) and the p.d.f. of k(n) is Pp(k(n)) = Pr(K = k) 
then by (4.74) and (4.75) 

P K (k(n)) = Pr[(f(n) - &)/h £ k(n) < (f(n) + l-j-)/h] 

= P r [hk(n) - |- < f (n) < hk(n) + |-] (4.77) 
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If the cumulative distribution function (c.d.f.) 
of f (n) is given by P p (f(n)) then 

P K (k(n)) = P p (a 2 ) - P p (a 1 ) ( ^ . 78 ) 



where 



a 2 = hk(n) - | (4.78a) 

a 2 = hk(n) + | (4.78b) 

This says that k(n) is distributed in such a way that its 
statistics are very similar to those of f(n) so that Equation 
(4.76) looks like it has four terms that are very similar 
in form to E[f (n-i ) f (n ) ] , the autocorrelation of f (n) . That 
is, it appears that the form of the autocorrelation of the 
error should be similar to that of the autocorrelation of 
the number to be quantized. 

In Appendix F it is shown that the autocorrela- 
tion of the state variables of a linear filter (i.e., an 
infinite precision linear filter) can be found, given the 
autocorrelation of the input and the transfer function between 
the input and the state under consideration. In the previous 
paragraph it is suggested that the autocorrelation of the 
errors is related very closely to the autocorrelation of the 
signal quantity to be quantized. In all cases this signal 
quantity is the product of a coefficient and either a state 
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or the input. The process of multiplying by a coefficient 
does not effect the autocorrelation of the signal quantity 
before or after the multiplier except to scale it by the 
square of the coefficient. Therefore, it is suggested that 
the autocorrelation of each error is similar in form to the 
autocorrelation of either a state variable or the input. 

A similar argument suggests that the cross correlation of 
error sources is related to the cross correlation of two 
signal quantities (or to the autocorrelation of a single 
signal quantity if the two coefficients multiply the same 
signal quantity). 

For simplicity, two extreme cases that bound the 
forms of correlation functions will be considered: (1) when 
the input is uncorrelated and (2) when the input is highly 
correlated (i.e., its autocorrelation is a constant). 

(1) Reduction for Uncorrelated Input. Appendix 
F shows that when the input is uncorrelated, the auto- 
correlation of the states will be a function bounded in 
magnitude by Ke“ Y ^ T where K is the variance of the state 
and y is proportional to the damping coefficient of the 
transfer function between the input and the state under 
consideration. But by Equations (3.9) &nd (3*11) it can be 
seen that the characteristic function for the state equations 
is the same as that of the overall transfer function so that 
y is the same for both states and output . 

The cross correlation between the uncorre- 
lated input and each state is proportional to the unit pulse 
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response of the respective state and the cross correlation 
between the two states is proportional to the correlation of 
the two unit pulse responses. All of these functions are 
then bounded in magnitude by an envelope proportional to 

e" Y M T . 

Assuming the correlations of the error 
sources behave in a manner proportional to the signals before 
the multipliers, these correlations may be written as 






K „e 
a£ 



— Y 1 i I T 



(4.79) 



where R a g(i) is defined as in (4.71f) and K a ^ is some positive 
proportionality constant involving the variance or covariance 
of the errors and ether factors . 

It may be possible to substitute (4.79) into 
a form of Equations (4.4l) which involves the magnitudes of 
the factors f^j and . Then to find |f^j| and |g^| one 
could examine 



<F> < <R (0)> + I {<A i ><R (i) > + [<A i ><R (i)>] t ) (4.80a) 

“ e i=l e e 



<G> < Z <A 1 ><R (i+1) > (4.80b) 

” i=0 ee 



<G’> < E <A 1 ><R ,(i)> (4.80c) 

” i=0 ee 



where < > indicates taking the absolute value of each entry 

of the matrix contained therein. This process may however 



give away too much to obtain a good estimate of the output 
error variance. This avenue is not pursued here but is 
suggested for possible future study. 

(2) Reduction for Highly Correlated Input. 
Appendix F shows that when the input is highly correlated, 
the autocorrelation of the states will also be constant. 
Assuming the error sources also have constant autocorrelation 

A A 

and cross correlation, the matrices R (i), R , (i) and 

R (i) in (4.21), (4.28) and (4.32) will be constant for all 
0 1 * 

i. For this case, let 



R e (i) = Q 



ee 



Q 



e l e l 




e 2 e l 



Q 



e l e 2 



Q 



e 2 e 2 



(4.8la) 



R ee - (i)=Q ee> = [Q e l£ ' 



R ee (i) = Q ee = «e 2 c 3 



( 4 . 8lb) 
(4 . 8lc) 



Then the matrix F whose entries appear in 
(4.48) can be written as 



00 oo . £ 

F = Q + Z k X Q a + Z [A X Q ] 
e i=l 6 i-1 6 

= Q e + A [I-A ] _1 Q e + { A [I- A ] _1 Q e } (4.82) 

°° n -1 

by the matrix identity Z A = [I-A] . Then 

n=0 
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tQ e,e (1+a ll" a 22“ a ll a 22 +a 12 a 21 ) + 2a 12 Q e 0 e, * 
f,. = — (1* . 83a) 



11 



1 + a + b 



tQe l e 2^ 1 a22 ^ +a 12 Q e 2 e 2 +a 21 Q e 1 e 1 + ^ a 22 _a ll a 22 +a 12 a 21^ Q e 2 e 1 } 



12 



1 + a + b 



(^4 . 83b ) 



tQ e 2 e 1 ^ 1 a ll )+a 12 Q e 2 e 2 +a 21 Q e 1 e 1 + ^ a ll a ll a 22 +a 12 a 21 )Q e 1 e 2 } 



21 



1 + a + b 



(i*.83c) 



{Q e 0 e 0 ^ a ll +a 22“ a ll a 22 +a 12 a 21^ + 2a 21 Q e,e * 
f 00 - 2-J. LJ_(i|.63d) 

1 + a + b 



But F must be a symmetric matrix, which 

implies that f 12 = f ?1 so that Equations (4.83b and c) imply 

that Q = Q making Q „ a symmetric matrix. This is 

e l e 2 e 2 e l 66 

consistent with the assumption that cross correlations are 

constant for all i, since when e contains R^Ci), Q e ^ e _^ 

contains R c (i) = R c (-i) = R ft (i). Then the entries of F 
(3a ot p up 

become 



{ ( 1+a, , -a^. o~a, a 05 +a, ) Q 



+ 2a, n Q 



11" 22 11 22 12 21 e 1 e 1 12 e i e 2 



} 



11 



-(4.84a) 



1 + a + b 



f = f 
12 21 



{a 21 Q e. ) e 1 +(1_a ll a 22 +a 12 a 21 )Q e 1 e 2 +a 12 Q e 2 e 2 } 



1 + a + b 



(4.84b) 



186 



(4.84c) 



22 



_ <2a 21 Q e 1 e 2 +(1 - a ll +a 22- a ll a 22 +a 12 a 21 )Q e 2 e 2 1 



1 + a + b 



The entries of the vector G (or G') which appear in (4.48) 
can be written as 

oo 

G = Z A 1 Q = [I-A]" 1 Q (4.85) 

1=0 

Since the object is to apply (4.84) and 
(4.85) to the canonic realizations of Chapter III, it can 
be assumed the Q = 0 since for all realizations either 

0 2^ 

e 2 (n) or e(n) is zero. Then 



g l = (1 - a 22 )Q e 1 e /<1+a+b) 



/I. A \ 

(4 . ooa; 



s 2 = a 21 Q eiE /(l+a+b) 



(4.86b) 



When Equations (4.84) and (4.86) are applied 
to the canonic realizations, then for the forms 



SM,TM: f._ = ( l-a-b)Q /(1+a+b) (4.87a) 

ll e 1 e 1 

f 10 = Q e ^ /(1+a+b) (4.87b) 

Id 

f 22 = 0 (4.87c) 

gl = g 2 = Q e /( 1+a+b) (4 . 87d) 



187 



SN,TN: 



(4.88a) 



SN,TN: f 


= -Q e e (4.88a) 

e l e l 


C\J 
i — I 


= Q_ _ /(1+a+b) (4.88b) 

e l e i 


C\J 

C\J 


= 0 (4.88c) 


e l = 


0 (4 . 88d) 



g 2 = Q e £ / ( 1+a+b ) (4.88e) 

since e 2 (n) zero f° r the untransposed forms. For the forms 



SM t ,TM t j f n 


= { (1-a-b )Q 0 + 2Q 0 o }/ (1+a+b) (4.89a) 

x 1 X 2 


f 12 


» {-b<3 + ( 1-b )Q 0 +Q a . }/( 1+a+b) (4 . 89b) 

e l e l e l e 2 22 


f 22 


= {-2bQ a +(l+a-b)Q ^ }/ (1+a+b) (4.89c) 

e l e 2 2 2 


S 1 = 


g 2 = 0 (4.89d) 


SN t ,TN t : f 


= {- ( 1+a+b )Q + 2Q }/( 1+a+b) (4.90a) 

e l e l 12 


f 12 


= {-( 1+a+b )Q ^ + ( 1-b )Q +Q e }/(l+a+b) 

e l e l 12 22 

(4.90b) 


f 22 


= {-2(l+a+b)Q . +( 3+a-b )Q }/( 1+a+b) 

e l e 2 22 

(4 . 90c) 
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0 



(4.90d) 



S 1 s 2 



These expressions can now be substituted 
into (4. 48). This yields for 



(c+erQ 



SM,SN,TM,TN: = 



e l e l 



2 ( c+e )Q 



( 1+a+b ) ' 



(1+a+b) 



e e 2 

+ (4.91a) 



SM t ,TM t : 



Av 



Q ^ + 2Q + Q 

e l 1 e l e 2 e 2 e 2 

(1+a+b) 2 



(4.91b) 



SN t ,TN t : 



Q 



e 2 e 2 



Av 



(1+a+b)' 



(4.91c) 



where, for the T forms, c and e are replaced by c' and e’ 
respectively. These may be written in another form by 
recalling that 



Q_ . = E[e . (n-i)e. (n) ] = E[e, (n)e. (n) ] 

e j j k j k 



(4.92) 



for all i. Then for 



SK.SN.TM.TN: o/ v = E[{. l4 : a + b - 



(c+e)e.(n) 0 

h + e(n)} 2 ] 



(4.93a) 



SM t ,TM t : 



e 1 (r.) + e ? (n) 2 

. 2 = e T (— } ] 

'Av U 1+a+b J 



(4.93b) 



SH t ,TN t : 



a 2 = E[(e„(n)/(l+a+b)} 2 ] 
Av 2 



(4.93c) 
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When the filter is in a limit cycle condition, 

it is hypothesized that a deterministic mean-squared error 

analysis would suggest replacing Q by the mean-squared 

e j e k 

error (effective value). Since it is well known that the 
limit cycles of the quantization error sources are bounded 
by h/2 for round off, then for that quantization method the 
rms value of each source must be less than or equal to h/2. 
Thus the rms value of the composite errors e-^, e ^ and e 
would have to be such that 




y 



2 h 

1 X 






e l e 2 1 



y l y 2 X 



(4 . 9*Jb) 




y 



2 h 

2 X 



(4.94c) 



e l e — 



h 

y l v X 



(4.94d) 




(4 .94e) 
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where the bar indicates "mean" and and v are defined 

after equations ( ^ .55) and (4.65a) respectively. Replacing 

the Q e e in (4.91) by the right hand side of the appropriate 

J k 

equation from (4.94) results in hypothetical bounds on the 
rms value of the output error as given in Table 4-2. These 
are conjectural and a detailed analysis is suggested for 
future research. 



One example of limit cycles given by Jackson 
[10] has been examined using the formulae of Table 4-2. The 
rms value of the largest limit cycle in the example was 
4.24h. The formula in Table 4-2 yielded a maximum rms value 
of 4.44h. This example is offered not as proof but rather 
as motivation. 



D. CONCLUSIONS 

Tables 4-1 and 4-2 list formulae for the output error 
variance for two extreme cases, uncorrelated error sources 
and maximum limit cycle rms values. It should be pointed 
out here that if any of the multipliers in the realization 
(not just in the transfer function) are such that it does 
not contribute quantization errors then these formulae need 
to be modified (by reduction of or v) to take this into 
account. For this reason, the expressions here cannot be 
directly compared to those previously presented in Table 2-1. 
Another reason is that previous analyses have not included 
any dependence on zero locations. This dependence is expli- 
cit in the expressions for the untransposed configurations 
since operations involving zeroes tend to follow those of 
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ARRAY 


M 1 


y 2 


V 


(Av) 2 

max 


SM 00’ SN 00’™0 TN 0 


2 


0 


2 


r l+a+b+c+e -i2 h 2 

L 1+a+b J “IT 


SM 01> SM 10- SN 01' SN 10'™1 


2 


0 


3 


r 2(c+e) + 3(l+a+b) -,2 h 2 
L 1+a+b J X 


SM 02’ SM 20 >SN 02 SN 20’™2 


3 


0 


2 


r 3( c+e ) + 2 ( 1+a+b ) -|2 h 2 
L 1+a+b J X 


SM 12 SM 21 SN 12 SN 21 


3 


0 


3 


A r l+a+b+c+e -| 2 h 2 
9C l+a+b ] X 


SM 11’ SN 11 


2 


0 


lj 


r 2(c+e)+^ (1+a+b) n 2 
L 1+a+b J X 


SM 22> Sfl 22 




0 


2 


r ^ ( c+e ) + 2 ( 1+a+b ) -| 2 h_ 
C 1+a+b ] T 


SM ij ™1J 


2 


2 


0 


16 [ 1 -] 2 h 2 

XDL 1+a+b J ip 


SN « 


2 


2 


0 


ur 1 i 2 h 2 
qL 1+a+b J X 



TABLE 4-2 MAXIMUM MEAN-SQUARED OUTPUT 
?RROR FOR HYPOTHESIZED 
LIMIT CYCLE ANALYSIS 
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poles in the flow between input and output in these realiza- 
tions, and thus are able to have a filtering influence on 
the internal noise generation. In the transpose configura- 
tions zero operations precede pole operations so the zero 
locations only influence the output error by the amount of ; 
correlation involved between error sources. 

An interesting observation regarding the expressions 
for error bounds by the several investigators is the recur- 
rence of various factors. For example, the factor, A, in 
Table 4-1 can be written 

A = (l-b)[(l+b) 2 - a 2 ] (4.95a) 

= (1-b) (1-a+b) (l+a+b) (4.95b) 

= ( 1-b ) ( 1- | a | +b ) ( 1+ | a | +b ) (4.95c) 

The factor (1-b) appears in the bounds of Sandberg and 
Kaiser, Long and Trick, and Yakowitz and Parker and in the 
estimate of Jackson. The factor (l-|a|+b) appears in every 
bound of Table 2-1 for real poles and in the single bound 
of Jackson and also of Bonzanigo. The factor A appears as an 
intermediate result in Parker and Hess [523 and Aggarwal [48]. 

Equations (4.10) and (4.21) are the most general result 
of this chapter. They could be very useful as a tool for 
understanding the interrelation among pole-zero location, 
structure and correlation in the accumulation of quantization 
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errors in the realizations of digital filters . Equations 
(4.48) are the general result for the realizations of Chapter 
HI* It is suggested that similar results could be found 
for the ladder structures by formulating their algorithms 
via state space techniques. Another area of promise is in 
that of multidimensional filters. State space formulation 
of these algorithms could possibly be expressed in terms of 
tensor analysis and the expressions of (4.10) and (4.21) 
likewise . 

Another area of application for these results is in the 
examination of realizations which perform quantization of 
signal variables after additions. In this case it is 
suspected that accumulated roundoff error will be reduced 
since the maximum composite error at time n would be one 
quantization step for truncation or half a step for roundoff 
rather than the sum of the maximum error from several sources. 
Quantization could also be performed before multiplications 
allowing double precision signal variables to flow as far as 
possible through the filter delays. There is often at least 
one path from the input to the output which w ould not require 
quantization . 

The results thus far still do not Indicate which is the 
"best 1 ' filter configuration. The reason is that the nature 
of the correlation between error sources has not been completely 
examined. Extension of the work of Appendix B is an avenue of 
endeavor which would be useful to pursue. Once this has been 
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accomplished, the "best" configuration for a given set of 
coefficients could be determined when the nature and 
statistics of the input signal are known. 

The results of this chapter can also be used in the 
determination of the optimal pairing of poles and zeroes in 
a cascade realization of higher order filters. 
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Figure 4-la. Realization of TMqq including quantization error sources 




Figure 4-lb. Realization of TM 0 J including composite error sources 
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Figure 4-2a. Non-Linearity of Rounding 
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V. DISCRETE FOURIER TRANSFORM FILTER 



A. INTRODUCTION 

In this chapter some relationships between the frequency 
domain and time domain representations of signals in finite 
impulse response (FIR) filters will be examined. These 
relationships provide insight and motivation for developing 
a new convolution generator algorithm which yields the cir- 
cular convolution of a given sequence (such as a finite 
impulse response) with sections of data taken from a continual 
data stream. Similarly, a new algorithm for performing the 
Discrete Fourier Transform of such sections of data is pre- 
sented . 
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B. DEFINITION OF TERMS 

Recall from section IIA4 that the circular convolution 

of two sequences of length N was defined by the sequence 

{y }, n = 0,1,...,N-1 such that 
n 



N-l 

y n ^ n u m h <n-m> 
m=0 



N-l 
Z u 
m=0 



<n-m> 



h 

m 



(5.1) 



where {u^}, n = 0,1,..., N-l are the two sequences to be 

convolved. The notation <•> indicates modulo N arithmetic." 1 ' 
Subscripts (or superscripts) are used as spatial or positional 
identifiers, whereas indices in parentheses will indicate 
time of occurrence or functional dependence. 

It was also shown that the same result (equation (5.1)) 
could be achieved via the equivalent operation of multiplica- 
tion in the frequency domain by defining: 

N-l 

U = E u W nK = DFT[ { u }] k = 0,1,..., N-l (5.2) 

k n=0 n n 

N-l ^ 

H = E h W n = DFT[{k }] k = 0,1,.. . ,N-1 (5.3) 

K n-0 n n 

Y k = U k‘ H k k = (5.*0 



Y To find <x> = x modulo y, write x = r + py such that 
P = [x/y] where [•] indicates the largest integer less than 
or equal to ( • ) . Then 0 < r < y and 
<x> = x modulo y = r = x - py . 
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where W = exp[-j 2tt/N] j (j = yj -l ) 

— 1 fh 

Note that W is the principal N root of 1. 



(5.6) 



The subscripts, k, above indicate the harmonic of 
interest, where the fundamental frequency is 



and T is the sample period in the time domain. 

In equation (5.5), if n is allowed to be outside the 
interval [0,N-1] then y^ = y <n >> i.e., the infinite sequence 



(y } is periodic in n with period N. 

CO 

C. THE DFT FILTER 

The term DFT FILTER, as used here, refers to a form of 
a finite impulse response (FIR) filter presented by Gold and 
Jordan [53]. This form will be examined in detail to gain 
knowledge about the nature of the output, thus providing a 
stepping stone to the development of a new convolution 
generator algorithm. 

1. Formulation 



filter, which is zero outside the interval 0 <_ n <_ N-l. This 
function has a Z-transform given by 



n = 2 tt/NT 



(5.7) 



CO 



Consider the impulse response, h(n), of an FIR 



N-l 

H(z) = I h(n)z 
n=0 



(5.8) 
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which is the transfer function of the FIR filter. Define the 
sequence, (h n > n = 0,1,..., N-l, such that h^ = h(n). 2 Then 
h(n) may be replaced by the Inverse Discrete Fourier Transform 
of its Discrete Fourier Transform, {H^} , as defined earlier. 
Then 



N-l N-l N-l N-l 

HU> ' n=0^ k-0 HkW >Z ' n = ^ E [H - Z (WVV} (5.9) 



N k^lf ‘ k n=0 



N-l , N 

Using the identity E x = — X | yields 

n=0 



H(z) = 1 " z 



1 - x 
-N N-l 



H 



k 



N 



ZT7 Ti 

k=0 1 - W z 



( 5 . 10 ) 



-IrN 

since W = exp[j2frk] = 1. 

Equation (p.IO) expresses the transfer function of the 
FIR filter in terms of the DFT coefficients, {H^}, of the 
sequence {h^} = (h(n)}. This is the form utilized by Gold 
and Jordan to show that "finite impulse response" and "non- 
recursive" are not synonymous since, as shown in Figure 5-1, 
each term of the summation in (5-10) is a recursive filter; 
yet H(z) represents a finite impulse response filter. 

But there is more involved in this form than is 
immediately obvious. 



O 

The dual notation is used in an effort to maintain 

consistency in the subscript/parenthesis notation. (h(n)} 

and {h } play a dual role here since (h } indicates coefficient 
n 

position in the filter, while, as an impulse response sequence, 
{ h ( n ) } is temporal in nature, yet they are identical in an 
FIR filter. 
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First, consider the terms, H^, in (5.10). Not only 
t h 

is H k the k harmonic of the DFT of {h(n)> but it is also 

the value of the transfer function, H(z), evaluated at the 
t h t h 

k N -root of 1, i.e., from (5.3) and (5.9) it can be 
seen that 



H k = H(z) 



z=W 



_ k - H(eJ k!!T ) 



NT 



(5.11) 



,-k . 



Second, note that, since W is, in general, a 
complex coefficient, the terms of the summation in (5.10) 
represent first order complex filters. Appendix shows 
that every first order complex-coefficient filter is also a 
second order real-coefficient filter. Furthermore, since 
|W" k | = 1, the poles of these filters lie on the unit circle 
in the z-plane; that is, each term represents a digital 
resonator with frequency 



_ , _ 2irk 

n k = kn = w 



( 5 . 12 ) 



(Cf. equation (5.7)). In other words, the form of the section 
of the filter in 5-1 between point 0 and the cross-section 
©-© is that of an oscillator bank and l-z N is a comb 
filter. The algorithm represented by equation (5.10) will 
henceforth be referred to as the DFT FILTER. The motivation 
for this phrase will become more obvious in succeeding sections 
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2 . Response to a Finite Duration Input 

It is well known that the response of any filter to 
an arbitrary signal is the linear convolution of that signal 
with the impulse response of the filter. The relationship 
between the response of an FIR filter to a finite duration 
signal (of length N) has an interesting connection to circular 
convolution. A more complete connection is made if one exam- 
ines the response of the FIR filter to a signal comprised 
of two successive copies of this signal of length N (total 
length 2N). Then a method for finding the circular convolu- 
tion can be devised. 

The response of any filter, H(z), to any input signal, 
u(n) is given by the linear convolution: 



When H ( z ) represents an FIR filter with an impulse 
response of length N, (h(n)j n = 0 , 1 , . . . ,N-1), and u(n) is zero 
except for n = 0,1,..., N-l, then by dropping all terms in 
(5.13) known to be zero (see Figure 5-2), 



oo 



y(n) = 2 u(m)h(n-m) 

m=-°° 



(5.13) 



r 



n 

I u(m)h(n-m) 
m=0 



0 < n < N-l 



N-l 

I u(m)h(n-m) N-l £ n £ 2N-2 



(5.1*0 



y(n) 




m=n-N+l 



0 



otherwise 
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since only on the heavy black portions cf the abscissa in 

Figures 5-2 d and e are the products in (5.13) non-zero. The 

only time (5.1*0 yields a term of the circular convolution 

of the two sequences (h(m)} = (h } and (u(m)} = (u } 

m m 

m = 0,1,..., N-l, is when n = N-l, since this is the only time 
when there are N products of terms taken from (u(m)} and (h(m)} 
in the summation. Then 

N-l 

y (N-l) = Z u(m)h(N-l-m) (5.15) 

m=0 



That this is a circular convolution term is assured by the 
fact that, in (5.1*0, because of the limits on m, (n-m) = <n-m> . 
But this is only one term of the circular convolution. 

To find the other terms, it is necessary to resort to another 
technique. By passing two copies of the data through the FIR 
filter all terms of the circular convolution are obtained. 

Since circular convolution can be viewed as linear 
convolution if one sequence is replaced by its periodic 
extension, consider the response of the FIR filter to an 
input comprised of two successive copies of the sequence 
(u(m)}, m = 0,1,...,N-1, i.e. 



r u(n) 

u(n) =, s u(n-N) 



0 < n < N-l 
N < n < 2N-1 



0 



otherwise . 



(5.16) 
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The response to u(n) would be (see Figure 5-3): 



C 



n 



Z •u(m)h(n-m) 
m=0 



0 < n < N-l 



(n) =< 



n 

Z u(m)h(n-m) N < n < 2N-1 
m=n-N+l 

2N-1 

Z u(m)h(n-m) 2N < n < 3N-2 

m=n-2N+l 



(5.17) 






otherwise . 



But in the second summation the limits on m imply 
that 0 <_ n-m N-l so it can be written immediately that 

n 

y(n) = Z u(m)h(<n-m>) N <_ n £ 2N-1 (5.18) 

m=n-N+l 

Furthermore, because of the constraints on n, the 
lower limit on m is always less than or equal to N, so (5.18) 
can be broken into two parts, yielding 

N-l n 

y(n) = Z u(m)h(<n-m>) + Z u(m)h(<n-m>) (5.19) 

m=n-N+l m=N 

N < n < 2N-1 



But since u(m) = u(m-N) for N <_ m <_ 2N-1 (from equation 
(5.16)), it follows that 
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N-l 

y(n) = E u(m)h(<n-m>) 

m=n-N+l 

N-l 

= E u(m)h(<n-m>) 

m=n-N+l 

N-l 

- E u(m)h(<n-m>) 
m=0 



n 

+ E u(m-N)h(<n-m>) 
m=N 

n-N 

+ E u(m)h(<n-m-N>) 
m=0 



N < n <_ 2N-1 (5.20) 



since <n-m-N> = <n-m>. The N outputs, y(n), occurring for 
N <_ n <_ 2N-1 are exactly the N members of the set of terms 
resulting from the circular convolution of the sequences 
(u(m)} and (h(m)} m = 0,1,..., N-l. Thus the fact that circu- 

j / 

lar convolution can be viewed as a linear convolution if one 
sequence is replaced by its periodic extension has been 
shown from a new perspective involving the FIR filter. 

But when one is processing continually received 
data, it is not efficient to stop and form periodic exten- 
sions of blocks of data to achieve circular convolutions on 
that block. 

Instead, the Fast Fourier Transform (FFT) has been 
used to perform high-speed circular convolution on blocks 
of data [ 15 - 19]. When these blocks are non-overlapping, 
this is a highly efficient method and there is little need 
to overlap when the signal being examined is highly stationary 
in a statistical sense. However, few real processes are 
highly stationary. To overcome this, overlapping blocks of 
data are processed. The question then arises, "How much 
overlap does one need and what will it cost in computation time?" 
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The following development provides an algorithm which yields 
the results in 2N-1 operations for complete overlap of data. 
This development for obtaining the circular convolution of 
blocks of data with the impulse response of an FIR filter is 
suggested by the symmetry property of convolution and the 
concept of periodic extensions of sequences. In equation 
(5.16) the periodic extension of (u(m)} was defined in order 
to find the circular convolution given by (5.20). Why not 
use the periodic extension of (h(m)}? Then the circular 
convolution of (h(m)} with (u(m)} could be found by the 
linear convolution of one copy of (u(m)} with the periodic 
extension of (h(m)}. Some very straight forward methods of 
doing this are possible. But with additional insight, a 
method of accomplishing this with the fewest number of summers 
and a simpler topology is possible. This is the object of 
the next two sections. The straightforward methods will be 
discussed after this development for comparison. 

3 . State Analysis 

The previous section states the objective of finding 
an algorithm which provides the circular convolution of a 
sequence (h(m)} with one block of data from a signal u(n) . 

It further suggests that some modification of an FIR filter 
would be suitable to this purpose. However, in that section 
only finite length signals were considered. From this point 
on, the objective is generalized to that of finding the cir- 
cular convolution of (h(m)} with every block of data of 
length N taken from an arbitrary continual data stream, 
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u(n) . Examining the states^ of the DFT FILTER will provide 
valuable insight into the nature of the signal quantities 
involved . 



Referring to Figure 5-1 
oscillator section as x k (n), k 
be the signal quantity at point 



, define the states of the 
= 0,1,..., N-l and let x(n) 




x(n) = u(n) - u(n-N) 



( 5 . 21 ) 



where u(n) is the input signal, which is assumed to be zero 
for n < 0 . 

For bookkeeping purposes, define the sequence 

{u (n)}, p = 0,1,..., N-l, (where n is now a time index), 

P / 

such that the elements of the sequence are 

u (n) = u(n+p-N+l) p = 0,1,..., N-l (5.22) 

These sequence elements may be viewed as samples 

taken from a set of functions u p (n) (see Figure 5 -4a for 

N = 4) also defined by equation (5.22). For example consider 

the sequence, {u (n')} = (3, 2, 1, 2} shown in Figure 5-^b 

P 

plotted versus p. The values 3> 2, 1 and 2 are the func- 
tional values of Uq(u), u^ (n) , u 2 (n) and u^(n) taken at 
time n' in Figure 5-4a. Notice that Figure 5-^b is identical 



^States are taken to be the signal quantites at the inputs 
to delay elements of the resonators of Figure 5-1. 
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to the region of the graph of u(n) (the top graph of Figure 
5-4a) indicated by a heavy black line on the abscissa. This 
heavy black line is the "window" for viewing the block of 
data referred to as (u (n')). The time index n* marks the 
right hand edge of the window. The left hand edge of the 
window is then n'-N+l. 

Next, define the sequence (s^Cn)}, k = 0,1,..., N-l 

as the DFT of (u (n)}. Note that the time index, n, associates 

P 

these two sequences in that (s^Cn)} is the DFT of the block 

of data in the window at time n, which is by definition 
h 

{u (n)}. Specifically, 

P . 

N-l .N-l . 

s, (n) = E u (n)W pK = Z u(n+p-N+l)W pK (5.23) 

K p=0 p p=0 

The relationship between the states of the DFT FILTER, 
x^(n), and the sequence (s^(n)} can now be found. 

From Figure 5-1, it is seen that 

x k (n) = VT k x k (n-l) + x(n) (5.24) 

By applying (5.24) iteratively (m-1) times it is 
seen that 

, m-1 

x. (n) = W kl V( n-m) + 2 W pk x(n-p) (5.25) 

* K p=0 



^It can also be seen that (u (n)> is exactly the set of 
quantities stored in the N-delay of Figure 5-1 a t tine n. 
u (n) is the quantity residing closest to but not yet arrived 
0 at the summer. The quantity at the summer is u(n-N) - u Q (n-i; 
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Assuming zero initial conditions, i.e. x.(n) = 0 
for n <_ -1, then for m = n+1, (5.25) becomes 



n _ , 

x k ( n ) = £ W-P x(n-p) (5.2 6) 

p=0 



Substituting (5.21) into (5.26) and assuming n >_ N-l, 



x. (n) = E u(n-p) W“ pk - E u(n-p-N) W~ pk (5.27a) 

K p=0 p=0 

N-l , n , n , 

= E u(n-p) W“ p + E u(n-p) W~ p - E u(n-p-N) W~ P 
p=0 p=N p=0 

Performing a change of variables in each of the three summa- 
tions of (5.27a), i.e. q = N-p-rl, q = p-N and q = p respectively, 

x (n) = ^"Vitn+q-N+l) W (q ~ N+1)k + E u(n-q-N) W“ (q+N)k 
K q=0 q=0 

- E u(n-q-N) W~ qk (5-27b) 

q=0 



Combining the second and third summations of (5.27b), realizing 
that W" (q+N)k = W" qk , yields 



x k (n) 



W 



N-l 

E u(n+p-N+l) 

p=0 




n 



E u(n-p-N) W pk 
p=n-N+l 

(5.27c) 



The summation in the first term of (5.27c) is, by 
(5.23), s. (n) . The second summation is zero by the assumption 
that u(m) = 0 for m < 0. So 
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( 5 . 28 ) 



x k (n) = W k s k (n) 

Equation (5.28) shows that for any n >_ N-l, the k" 

state of the DFT FILTER is proportional to the k th harmonic 

of the DFT of the sequence, (u (n)}, that is, of the sequence 

P 

consisting of the present input and the past N-l inputs. (In 
fact, this is true for n < N-l if one considers sequences 
with leading zeroes.) Thus (5.28) provides the motivation 
behind the name DFT FILTER for the algorithm of equation (5.10) 
and Figure 5-1. 

The state response of the DFT FILTER for two specific 

types of input sequence will now be considered as examples 

to illustrate the operation of the comb filter section. 

a. Example 1! Unit Pulse Response 

Let u(n) consist only of the unit pulse at n = 0. 

This pulse is sent directly to the resonators as x(0) and 

starts their oscillations. It is also stored in the N-delay 

of Figure 5—1. (Figure 5—5 shows the result of this input 

— 1c n 

for N = 4.) These oscillations are such that x R (n) = (W ) 

n = 0,1,... ,N-1 . 

From (5.22), u Q (N-l) = u(0) = 1 and u (N-l) = 0 
for p = 1,2,..., N-l. So { Up (N-l) } is an impulse sequence 
with the impulse at p = 0. On the other hand, x k (N-l) is 
equal to (vr k ) N-1 = W k as expected from (5-28) since the DFT 
of an impulse sequence is {s^} = (1^ • The factor V' can be 
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viewed as a phase delay of 2Trk/N since W k = exp [-j 2TTk/N] . 

This phase delay is evident in Figure 5-5. 

Recall that the original pulse was stored in the 
N-delay and now impinges on the resonators as a negative 
pulse at time n = N since x(N) = u(N) - u(0) = -1. But 
x^(N) = W k x^(N-l) + x(N) = 0 and so for n _> N the states 
of the DFT filter will be zero. 

In effect then a unit pulse input causes the 
resonators of the oscillator section to oscillate for one 
period, n = 0,1,..., N-l, after which time they are turned 
off. Since an arbitrary input is simply a string of weighted 
and delayed impulses and since linearity applies, the state 
response of the DFT FILTER to an arbitrary input is the 
superposition of a string of weighted and delayed oscillation 
periods . 

Before proceeding with example 2 some additional 
insight into the intimate relationship between the states 
of the DFT FILTER and the DFT of sequences will be helpful. 

Consider equation (5 *23) which by a change in variable 
(p is replaced by p-1) can be rewritten as 

sAn) = lu n (n) W (p_1)k (5.29) 

k p=l p_i 

From (5.22) it follows, in general, that 
u (n) = u( (n-m) + (p+m)-N+l) = u (n-m) (5-30) 

p . 
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for any m. So (5.29) becomes 



s (n) = W _k ® Up(n-D W Pk 
P = 1 
. N-l 

= W" K { Z u (n-l) W pk + u M (n-l) 

p = 0 P N 

= w" k (s k (n-l) + u N _ 1 (n) - u 0 (n- 

- W k (s k (n-l) + u(n) - u(n-N)} 

= W k (s k (n-l) + x(n) } 

Substituting (5.31) into (5.28) yields 

x k (n) = s k (n-l) + x(n) 

Equation (5.32) says that x k (n) is also 
the DFT of {u p (n-l)} plus an additional term. 
(5.31) iteratively (m-1) times it can be shown 

m 1 

s, (n) = W"* s.(n-m) + Z v/~^ q+1 ^ k x(n-q) 

K k q=0 

Substituting (5.33) into (5.28) yields 

x (n) = w -(m-l)k s ( n _ m ) + z W _qk x(n-q) 

K K q=0 



W Nk - UQ(n-l)W^) 
1 )} 



(5.31) 



(5.32) 

equal to 
By applying 
that 

(5.33) 



(5.3*0 
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So the states of the DPT filter at time n contain a phase 
shifted version of the DFT of previous sets of window data 
plus some additional terms. Example 2 will demonstrate 
this with a periodic input signal. 

b. Example 2 ' Response to Periodic Input 

Let u(n) be zero for n < 0 and periodic with 
period N for n _> 0 (See Figure 5-6). Then x(n) (point 0 
in Figure 5-1) consists of only one period of u(n) while 
the window of data stored in the N-delay after n = N-l 
looks like a continually progressing circular shift of that 
period. The DFT's of such circularly shifted sequences 
is known to be identical except for a phase shift. Since 
x(n) = 0 for n _> N the resonators will continue to oscillate 

rr-|V> qyic "P r \ yr \ q +“ o v-i +- mo n f 1 1 r) o Vnif nKoMfri nrp nkoop 

w AiUiliO ■w 1 ^ J. U U-i i 1/ ImwV jJ. . -L. 04. V4. v-t W/ 1-4, i. 1 — 1 1 

at the rate, 2irk/NT radians/seccnd . This is predicted by 
equation (5.3*0 when m is not so large that x(n-q) is a 
term from the single period of data which enters the 
oscillator portion of the DFT FILTER. 

In the next section significant use will be made of 
the fact that x, (n) contains some information about the 

ri 

DFT’s of any number of sequences preceeding (u (n)} in time. 

Ir 

That this information is contaminated by phase shifts and 



additional terms will not prove insurmountable. 

^ • Output Analysis 

The relationship between the states of the DFT FILTER 
and the DFT o'f the input sequences provides, as shown below, 
the knowledge which makes possible the extraction of complete 
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information about the circular convolution of every block 
of data of length N with the impulse response of the FIR 
filter. 

According to Figure 5-1, 

! N-l 

y(n) =jr EH x (n) (5-35) 

N k=0 K K 



But in the previous section (equation (5.28) it was 
shown that x^(n) was proportional to s^(n) . Thus, 



y(n) - i Vh k W k . k (n) = J JV . k („) W< N - 1)k 

k-° k=0 <5.36) 



which is, from equations (5.5) and (5.1) the (N-l)th term 

of the circular convolution of the sequence {h } with the 

n 

sequence {u (n)}, since (s, (n)} is the DFT of {u (n)}. 

P k p 

So the output may be written as 



N-l N-l 

y(n) = l u <N-l-m>^ = 2 u m^ n ^ h ( <N-1_m> ) (5.37a,b) 

m=0 ~ m=0 



or 



N-l N-l 

y (n) = l u M , (n) h(m) = I u_(n) h(N-l-m) (5.37c,d) 

m=Q N " i-m m=0 m 



since <N-l-m> - N-l-m for the values of m allowed in the 
summations. These four expressions are equivalent. 
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However, if equation (5.32) is substituted into (5.35), 

then 



. N-l N-l 

y(n) = i E H s (n-l) + ± Z H x(n) (5.38) 

w k=0 K K N k=0 k 



The first summation is the zero term of the circular convolu- 
tion of {h^} with {u (n-l)}. Since x(n) does not depend on 
k, the second term of (5*38) is merely x(n).h(0) since 

is, by definition, the zeroeth term of the IDPT of 
{H^}. Therefore, not only is y(n) the (N-l)th term of the 

{u (n)} convolution, it is also possible to extract the 
P 

zeroeth term of the convolution with (Up(n-l)} by subtracting 
x(n) h(0) from y(n). In fact, N-l more convolution terms 
can be extracted from y(n), as will be shown next. 
Substituting (5-3^) into (5-35) yields 



y(n) 



N-l 

= M ^ H k s (n-m) 
N k=0 K K 

, N-l 

= N 2 H k s k (n “ m) 

N k =0 K K 



w -(m-l)k 

w -(m-l)k 



+ 



1 

N 



N-l 
Z H 
k=0 



m-1 
Z W 



-qk 



q=0 



m-1 

+ Z x(n-q) h (q ) 
q=0 



x(n-q) 

(5.39) 



the first term being the (m-l)th term of the (Up(n-m)} 

convolution. For notational purposes define y^(n;k) as the 

mth term of the circular convolution of the sequence 

\ 

(u (k) } with (h }, where the index n indicates when that 
P P 

term is available from y(n). V/ith this notation then, 
equation (5.39) says 
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(5.40) 



m-1 

y -.(njn-m) = y(n) - Z x(n-q) h(q) 

q=0 

This algorithm is shown in Figure 5-7, using the DFT 
FILTER form of H(z). Figure 5-8 uses the FIR canonic form 
of H(z) and utilizes the normal delays of that form in 
conjunction with a direct feed-forward path to form 
x(n) = u(n) - u(n-N) . 

The algorithm of equation (5.40) might be called a 
"circular convolution generator." It is most efficient when 
one must calculate all N convolution terms for every subse- 

oo 

quence of length N of the "infinite" sequence {u(n)} Q . In 
that case, it requires, on the average, roughly two multi- 
plications per term (i.e., 2N multiplications per convolution) 
instead of the N multiplications per term required to perform 
each convolution separately or the (log 2 N) + 1 multiplica- 
tions per term required to accomplish each convolution via 
the FFT. 

It may be noted here that it is not necessary to 
generate the term y^ -^(n;n-N) since this term is available 
as y(n-N) = y^ 1 (n-N;n-N). This saves one delay, one 
multiplication and two additions but is insignificant for 
large N. 

To illustrate when the terms of the convolution of 
a particular sequence (Up(rn)} become available at the 
outputs of the circular convolution generator of Figure 5-7 
or 5-8, Figure 5-9 shows the time base of the N outputs. 
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The points connected by the diagonal line represent the N 

terms of the circular convolution of {u (m)}. The last 

P 

output is the redundant term y N (n;n-N). 

Section B alluded to some straight-forward methods 

of performing this same operation. Figures 5-10 and 5-11 

are two such methods. The sequence (u (m)} is stored in 

P 

the first N delays of these algorithms at time n = m, there- 
fore these delay elements are equivalent to the N-delay 
window positions of the DFT FILTER. So by the notation 

■f- 

convention, y m (n;k) represents the nr term of the circular 
convolution of the data in the first N delays at time k and 
is available at time n. 

Figure 5-10 requires 2N-2 delays, 2N-1 multipliers, 

and N(N-l) summers not including the redundant operations 

at the right end. The worst thing about Figure 5-10 is 

the horrifying number of circuit crossings required. 

Figure 5-11 reduces the number of summers and circuit- 

crossings by recognizing that to achieve each successive 

output term moving from left to right it is only necessary 

to add one product and subtract another. Thus Figure 5-H 

requires 2N-2 delays, 2N-1 multipliers and 3N-5 summers. 

Again the redundant operations have been removed from con- 

2 

sideration. But there are still N - 2N +3 circuit 
crossings ! 

Figure 5-12 is a rearrangement of Figure 5-8 which 

5 

has no circuit - crossings . This arrangement of the circular 
convolution generator requires 2N-2 delays , 2N-1 multipliers 
^This is particularly advantageous in LSI technology. 
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and 2N-2 summers. It is believed that this is the least 



number of operations required to obtain the objective. 
Obviously the topological circuit crossings cannot be 
reduced further. 

Note that the portion of Figure 5-12 to the left of 



cross-section 



©-0 



could be replaced by a DFT FILTER 



but to no apparent advantage at this time. 



D. MODIFIED DFT FILTER 

In Section C3, it was shown that the states of the DFT 
FILTER were proportional to the DFT of the input sequence, 

(u (n)}. It is also possible to design the oscillator 
section so that the states are exactly the DFT of {u (n)}. 

s' 

This will lead to a minor refinement of the convolution 
generator. It will also lead to some observations that will 
be pertinent in Section E where a new DFT algorithm is 
derived . , 

Taking the Z-transform of equations (5. 24) and (5.28), 



X k (z) = 



X( z) 



n H“ k - 1 

1 - W z 



= VT 



V z) 



k = 0,1,... ,N-1 

(5.41) 



whe re 



X k (z) = Z{x k (n)} 


(5.41a) 


X(z) = Z { x ( n ) } 


(5.41b) 


S k (z) = Z{s k (n)} 


(5.41c) 
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The right-hand equality of equation ( 5 . 4 1 ) may be 
rearranged to provide a transfer function between X(z) and 
S k ( z ) , i .e . 



S R (z) w -k 

nvr - x . w-v 1 



(5 .* 12 ) 



The right hand side of (5. *12) represents a digital resonator 

_ \r _ \r fs 

with a feed-forward gain of W and a loop gain of W z 
Figure 5-13 shows a modified DFT FILTER where (s^Cn)} are 
the states of the oscillator portion. 

Designating the output of the modified DFT FILTER as 
v(n), the transfer function of this form is 



rt f \ 



v(z) = i - 
uTz) ' n 



-N N-l W“ k H, 
z v k 

-k -~\ 

k=0 1 - W z x 



(5.*»3) 



In the form of an FIR transfer function, (5.^3) becomes 



N-l 

G(z) = I h(n+l)z (5.***0 

n=0 

where h(N) is by convention equal to h(0) to be consistent 
with the concept of periodic extensions of sequences. It 
would also be possible to replace h(n+l) by h( < n+l > ) which 
is consistent with the "modulo N" nature of circular convolution. 



6 This form is- related to the Goertzel algorithm [ *1 ] for 
finding the DFT (see Section HA 1 ! ) . 



220 



Upon examination of the output of this filter, it is 
seen that 



i N-l 

v(n) = w l H s (n) ( 5 . 115 ) 

k=0 K K 



which is the zeroeth term of the circular convolution of 
(Up(n)} with 

Substituting equation (5-33) into (5.^5) , 



T N_1 mlr i N-l m-1 / ,, x. 

v(n) = i EH, s (n-m) W" mk + i- E H, E W (q+1)k x(n-q) 
k=0 K w k=0 K q=0 



m-1 

= y (njn-m) + l x(n-q) h(q+l) 
m q=0 



(5.46) 



This shows again that the output contains sufficient informa- 
tion to allow the extraction of N convolution terms in a 
similar manner as was done previously. In this case 



m-1 

y (n;n-m) = v(n) - I x(n-q) h(q+l) (5.47) 

m q=0 

This algorithm is shown in Figure 5-14 with G(z) repre- 
sented in the FIR form. G(z) may also be implemented, of 
course, in the modified DFT FILTER form. In either case. 
Figure 5-14 may be referred to as the "modified circular 
convolution generator." 

The modified convolution generator has no particular 
advantage over the one derived in Section C except that the 
temporal order in which the terms of a convolution associated 



with a particular sequence {u (n)} become available is also 

P 

the ascending order of the index, m, in {y (n; k)} , 

m = 0,1,...,N-1. Furthermore, this result occurs without 

relying on the redundant multiplication noted in Section C. 

Note also that instead of the apparent circular shift in the 

index, m, there is, in Figure 5-6, a circular shift of the 

index on the coefficients, h , and that the redundant multi- 

n 

plication would be by a second h^ which does not appear. 

Extending this idea, it is then possible to extract the 

convolution terms in any circularly shifted order by a 

suitable circular shift of the h 's. It will be shown in 

n 

Section F that by relying on the symmetry property of con- 
volution it will be possible to extract the convolution term 
in any circularly shifted order by a suitable circular shift 
of the Up's. This will have a significant impact on a possible 
application of the circular convolution generator. 

E. A NEW DFT ALGORITHM 

In the previous sections, the primary interest was in 
a time domain extraction of circular convolution, though 
frequency domain techniques were resorted to in order to 
gain an insight into the properties of the available signal 
quantities. In this section, a new DFT algorithm is derived 
by making use of some of those same insights. 

In Section D, a transfer function between X(z) and S k (z) 
was derived which provided the spectrum of the sequence 
{Up (n) } at the output of the oscillator section (cross 
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section 



©-© 



of Figure (5-13)). This part of the 
algorithm is not new. It has been derived previously by 
Dillard [54] as a means of generating spectra 
or spectral components of a series of sequential sequences 
{u (n)}. as defined in Section C3. This algorithm may be 
described as a moving (sliding) window DFT or uncoupled DFT. 

As a moving window DFT, it may be observed that one may 
feed one sequence of length N into the resonators and obtain 

p 

the entire spectrum of that sequence in N multiplications, 
the same number required to perform the DFT in its definitive 
form. However, the algorithm requires only N multiplications 
per spectrum, on the average, if one is interested in full 
overlap of the moving window. Further savings are realized 
if one is interested in only a few spectral components, since 
the algorithm provides each harmonic independently of all 
others, i.e., it is an uncoupled form. 

This algorithm has one major disadvantage which has not 
been expounded or examined to date. Theoretically this 
algorithm provides the exact discrete spectrum of each 
sequence, (u p (n)}, independent of every other (u p (m)}. How- 
ever, when finite precision arithmetic is used, significant 
errors destroy that idealistic situation. 

First of all, it is apparent from studies of possible 
pole locations in second-order, finite word-length, real coef- 
ficient digital filters [11 ], that no method has been de- 

vised (or can be devised) to place the poles of a set o± such 
filters exactly on the unit circle in the Z-plane and evenly 



223 



spaced in frequency (except for the set {z|z = 1, j, -1, -j}) 
because of the finite length of the coefficient representation. 
This means that there will be errors due to the fact that the 
transfer function realised is not the one required. Prom 
another point of view, the problem is that W nk in the defini- 
tion of the DFT is not representable, in general, in a finite 
word length. So this difficulty applies to any realization 
of a DFT algorithm. 

Secondly, in a finite word length algorithm, it is also 
necessary to quantize the results of multiplication because 
of the extended word length that operation implies . This 
difficulty, too, applies to any realization of the DFT. 

But, the third and most significant problem to be faced 
here is the recursive nature of the comb filter section of 
Figure 5-13* In effect, it magnifies the problem of quanti- 
zation after multiplication. 

Recursive digital filter sections have been extensively 
studied in the literature [ 39> 51 ] and in Chapter IV of 

this work with respect to error generation and propagation. 

In particular, recursive filters are highly subject to the 
generation of limit cycle conditions, especially when the 
pole locations are close to the limits of stability . In 
the case of the DFT FILTER the pole locations are exactly 
on the unit circle, the stability boundary. It should suffice 
here to say that recursion should be avoided whenever possible 
or practicable . 
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There is a set of complex coefficients, however, which 
do not contribute errors of the type described above. These 
are the coefficients {1, j, -1, -j } . If a transfer function 
contains these and only these coefficients, there will be 
no shift in pole location between the required and the 
realized filters. Furthermore, there will be no error intro- 
duced by quantization after multiplication by these coeffi- 
cients and hence no recursive propagation of such errors, 
even in a recursive filter. This seems to be a trivial idea, 
but is importance comes to light in view of an observation 
by Hess [11] regarding the realization of digital resonators 
of a desired frequency. 

Hess observed that, since the frequency of a resonator 
is dependent upon the sample period, any frequency can be 
generated. By choosing to^T equal to an integer multiple 
of tt/ 2, only multiplications by 1, ], -1 or -j are required 
to have an oscillator with frequency as shown below. 

It was shown in Section D that the DFT of windowed data 
could be achieved through an appropriate arrangement of 
digital resonators whose poles were located at the associated 
harmonic frequencies . An algorithm will now be developed 
which exploits the relationship between frequency and sample 
period. 

Consider the complex resonator with the transfer function 

H.'U) = i— r (5. 'i8) 

1 1 - jz 
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H(z) has a pole at 



z = j = exp[jir/2] 



(5.-49) 



Equating the right hand side of (5.49) with the exponential 
form of z on the unit circle: 



z = exp[jco d T] 



(5.50) 



yields the resonator frequency 



IT 

W d 2T 



(5.51) 



Furthermore, one may also observe that this pole is at 



z = W 



N 

T 



(5.52) 



which means that this resonator would generate the DFT 
frequency 



o _ Nfi N 2tt 
W N/4 T T * NT 



(5.53) 



and this is exactly in equation (5.51). 

Now, if T were to be replaced by T' = 2T then this pole 
location would correspond to half the original frequency . 

However, it 'is not part of the approach at this point to 
consider changing the sampling rate, since it is assumed 
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that data is being received at a fixed sampling rate. On 
the other hand, consider the complex resonator 



H 2 (z) = 



1 - jz 



-2 



(5.5*0 



and evaluate this function on the unit circle. Thus 



jw.T 

H 2 (e d ) = 



-ju 2T -jtu.T' 

1 - je d 1 - je d 



(5.55) 



This is exactly H^z) evaluated on the unit circle when T 
is replaced by T' = 2T. The sampling rate had not been 
changed in equation (5.55) but the effect is similar. In 
fact, if H 2 (z) were used in a desampling mode (where only 
every other output is considered) the effect would be the 
same . 

But now consider the poles of H 2 (z) which occur at 



z = 




±[e j7T/2 ] 2 = e J7T/i) , 




(5.56) 



One could achieve the effect of a single pole at 
z = exp[jfi/*l] by modifying H 2 (z) to include a zero at 
z = exp[j 5^7*1 1 • Thus 



H 2 (z) 



i - z - 1 = i - vr^/ 8 z- 1 

1 - jz” 2 1 - jz -2 



(5.57) 



22 ? 



The pole, z exp[j7r/^] - W ^ 8 , Is one required in any 
DFT FILTER algorithm which has N equal to an Integer 
multiple of eight . 

Immediately obtainable from H 2 (z) is another pole required 
in multiple-of-eight algorithms; i.e., z = W“ 5N/8 . In fact 
all eight frequencies for N = 8 are obtainable from the set 
of transfer functions: 



G N/8 (z) 



G 5N/8 (z) 



G 2N/8^ 



G 6N/8 (z) 



G 3N/8 (z) 



! - W - 5 N / 8 Z “1 



—^2 = G ± (z) = H 2 (z) 



1 - jz 

1 - W“ N/8 z- 1 



1 - JZ 



-2 



- G(.(z) 



, ..-6N/8 -1 

1 - W z 



-2 



1 + z 



1 ■ 

1 + z" 2 

1 - W-7M/8 s- 1 

1 + jz" 2 



- G 2 (z) 



- Gg(z) 



- G^(z) 



G 7N/8 (z) 



G 8N/8 (z) 



G 4N/8 (z) 



- * " . G_ ( s , 

1 + jz 1 

, „-4N/8 -1 

_ 1 - W z _ P 

g q (z) 

1 - z 

= 1 ■- — = G^(z) 

1 - z 



(5.58a) 
(5.58b) 
(5.58c) 
(5.58d) 
(5.58e) 
( 5 . 5 8 f ) 
(5.58g) 
(5.58h) 



where the subscript indicates the harmonic number of the 
associated frequency. That is, the subscript corresponds 
to k in Figure 5-13. 
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Figure 5-16 depicts an algorithm which makes use of this 
development for N = 8. As in Figure 5-13 the input must 
first pass through the section (1 - z _N ) before reaching 
the resonators. Also, recalling that each resonator in 
Figure 5-13 had a feed forward gain of W~ k , each forward 
path in Figure 5-16 includes this scale factor which was 
not included in equations (5.58a-f). 

Figure 5-16 is a rather complicated algorithm for an 
N = 8 DFT. It also appears to require twice as many 
multiplications as the algorithm of Figure 5-13. Some 
further refinements would make this a more desirable form. 
Examination of equations (5.58) yields these refinements. 

Consider equation (5.5 8h ) . Multiplication by W° is 
trivial since W G - 1. In addition, the denominator of 
( 5 . 58 h) is factorable into (1 + z~ (1 - z ^) which means 
(5.58h) can be rewritten as 



G lj (z) 



1 + z 



-1 



(5.59a) 



In (5.58g) a similar result occurs. 



Here W 



-W8 



- 1 , 



so that 



G 0 (z) = i - z -1 



(5.59b) 



In (5.58c), V/ = j; thus 



r ( s 1 - jz -1 L 

c, 2 (z) ~ 7^ ; . _:r, , . _ - 1 ■ 



(1 + Jz" ) (1 - Jz ) 1 + jz 






(5.59c) 
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And in (5.58d), W" 2N ' /8 = -j ; so that 



C 6 ( z) 




(5.59d) 



Equations (5.59a-d) denote resonators with a single delay 
element, rather than the two delay elements in equations 
(5.58 c,d,g and h) . The cancellations which were used to 
derive equations (5.59a-d) from (5.58 c,d,g and h) are not 
possible in equations (5.58 a,b,e and f) without reintro- 
ducing non- trivial multiplications into the recursion loop. 

(A non-trivial multiplication is one where the coefficient 
is other than 1, j, -1 or - j . ) 

Another simplification involves flow diagram manipulations 
to be performed on the zeros associated with the resonators 
of equations (5.58 a,b,e and f) including the scale factor 
W “ associated with the appropriate harmonic. 

Figure 5-17a shows the zero paths associated with 

-5N/8 

equations (5.58 a and b) after observing that W is 

equal to -W”^ 8 . (In fact, this is true independent of the 
value of N, since 



W -^ N/8 = exp[(-j2TT/N)(~W8)] = expCjTt] = -1.) 

( 5 . 60 ) 

It is now possible to remove the multiplications by 
W N/8 which are closest to the output in Figure 5-I7a by 
placing a multiplication by V/ * ^ 8 in each input to the 
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summers as shown in Figure 5-17b. Then one can consolidate 
consecutive multiplications by W -1 ^ 8 to obtain a multipli- 

O XT / Q 

cation by W which is actually a multiplication by j 

as shown in Figure 5-17c. The next manipulation is to 
remove identical multiplications which occur in all paths 
leaving a node by placing a single multiplication before 
that node. Figure 5-17d shows the final arrangement of a 
flow chart for equations (5.58 a and b) which is equivalent 
to that shown in Figure 5-16. (Residual multiplications by 
-1 have been moved to the right and made part of the summer 
(subtracter) operation.) A similar manipulation is possible 
for the zero paths of equations (5-58 e and f ) . 

Figure 5-18 depicts an equivalent flow diagram to Figure 
5-16, incorporating all the above-mentioned modifications. 
(Note, only non-trivial multiplications are drawn inside a 
square. All trivial multiplications, i.e., multiplications 
by 1, j -1 or -j , are shown inside circles.) Three major 
properties of this figure should be observed. 

First, all non-trivial multiplications have been left 
in the form W” kN//8 , because in this form they are independent 
of N. On the other hand, the subscripts on s^(n) are also 
shown in this manner to indicate their dependence on N. The 
usefulness of this fact will be seen later. 

Second, Figure 5-9 is a somewhat uncoupled form. The 
four single-delay stages are completely uncoupled and the 
two double-delay stages are pairwise uncoupled. 
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Third, there are no non-trivial multiplications in the 
recursive part of the algorithm. Therefore, no recursive 
error propagation will ensue . Errors generated in this 
operation are only in the output zero-paths and only effect 
the DPT of one window of data. 

A final observation on this specific example is that the 
number of non-trivial multiplications involved is reduced 
from what was required in Figure 4. There are six nontrivial 
multiplications required in Figure 5-13 For N - 8, while 
only two exist in Figure 5-18. A first estimate on the 
number of multiplications for the FFT is N log 2 N = (8) (3) = 
2^, but through manipulations of flow charts similar to that 
performed on the zero paths; the FFT algorithm can also be 
reduced to two non-trivial multiplications. So nothing has 
been gained in terms of non-trivial multiplications over 
the FFT, but some degree of uncoupling has been achieved. 

This will shown to be an advantage in succeeding paragraphs. 

When N = 1 6 , four resonators (with four delays each) are 
required to provide the 16 poles necessary for the transfer 
functions , 



w -kN/l6 

G kN/l6 (z) ~ 1 _ w -kN/lb z -l 



(5.61) 



Again the feedback coefficients for the four resonators will 
be 1 , j , -1 and -j . 

It is necessary in this case to cancel three poles of 
each resonator by zeroes to provide the single pole of (5«6l). 
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In the process of doing so it can be observed that the 
resonators with feedback coefficients 1 and -1 reduce 
exactly to Figure 5-18 with N = 16 . (Notice that N = 16 
does not change the coefficients but only changes the sub- 
scripts on the outputs , since now W = exp[-j 2 tt/N] = exp[-j 2 tt/i 6 ] . ) 
This reduction occurs through the same reasoning by which 
Equations (5.59) were derived. 

The remaining two resonators and associated zero paths 
are shown in Figure 5-19. When the appropriate flow diagram 
manipulations are performed, Figure 5-19 reduces to Figure 
5-20. Figures 5-18 and 5-20 then, taken together form a 
complete DFT algorithm for N = l6 . The outputs of Figure 
5-l8 are the even harmonics and those of 5-20 the odd ones. 

If the algorithm of these two figures taken together v. r ere 
to be built or programmed, there would be available to the 
user DFT algorithms for N = 1, 2, *1, 8 or 16 . The only 
parameters to be chosen are: 

(1) the length of the N-delay at the front end 

(2) the number of resonator loops to take advantage of 

(3) the subscripts of the outputs 

But all of these parameters depend on N and can be incor- 
porated in one switch or program statement which indicates 
the value of N. Notice also that it is only necessary to 
store two non-trivial complex coefficients since 

w -6n/16 0 w - 3N/8 . jW -N/8 (5.62a) 
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(5.62b) 



W “ 9N/l6 = _ w - N/l6 



and 



w - 2N/l6 = w " N/8 



(5.63c) 




be stored. These can further be reduced to two real 
coefficients since 



Another reduction of the algorithm is possible if only 
real data is being processed. Then all the resonators with 
-j as the feedback coefficient can be removed since the 
outputs associated with them are available as the complex 
conjugates of the outputs of the resonator with the same 
number of delays and feedback coefficient, j; e.g. 

S 7N/I6 (n) = S 9N/I6 (n) (5 ' 65 

where * indicates complex conjugate. The output 3 yjj/^g( n ) 
is in Figure 5-20b while s^^g is in Figure 5-20a. 



R [W N/8 ] = -R [W“ 3N/8 ] 

c G 



(5.64a) 



i rw“ N/8 i = I [W 3N/8 1 

m* - m 



(5.64b) 
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Algorithms for N equal to a power of 2 are possible 
for all higher powers as well. Figure 5-21 shows an eight- 
delay stage with feedback coefficient, j , required for N = 32 
and higher powers of two. Figure 5-22 shows the equivalent 
stage after reduction of the number of non-trivial multipli- 
cations. Figures 5-23 and 5-24 show the sixteen delay stage 
required for N = 64 and higher powers of two. 

Only the stages with feedback coefficient j are shown 

since the form for the stage with -j is the same except that 

-N/32 

all powers of W J in Figure 5-22 are replaced by the same 
powers of v/”^N/32 and all powers of w“N/64 in Figure 5-24 
are replaced by the same powers of w“3N/64. 

Taking Figures 5-18, 5-20, 5-22, and 5-24 together as 
one algorithm, there are available the routines for finding 
the DFT with N = 1, 2, 4, 8, 16, 32, and 64, and the only 
parameter that needs to be changed is N. 

Another advantage of this algorithm is that the semi- 
uncoupled nature of it allows one to reduce the number of 
calculations required when only certain subsets of the har- 
monics are required. For example, suppose one is interested 
only in the odd harmonics of real data for N = 64 , Then 
only the algorithm of Figure 5-24 need be performed (inclu- 
ding .ie feedforward (1-z path of Figure 5-18) . If only 
every other even harmonic beginning with k = 2 for N - 64 
is desired, then Figure 5-22 yields all the required 
information, when (1-z ) precedes it. 
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If the data being processed are not from a continual 

data stream, but rather are just N pieces of data, then the 

-N 

delay z is not necessary. The data can be fed sequentially 
to the resonators as long as the delay elements have been 
cleared (set to zero initial conditions) before entering 
the first datum. The calculations on the zero paths need 
only be performed at time n = N (if it is assumed that the 
first datum is entered at time n = 1). 

From the previous paragraph and examination of the 
resonators of Figures 5-18 through 5-24, it becomes clear 
that this method of finding the DFT is essentially different 
from other methods since the first operations performed on 
the data involve adding all data, every other datum, every 
fourth datum, etc., after appropriate multiplication by 
powers of 1, j, -1, -j . That is, if this method were performed 
in parallel from the start rather than with a sequential data 
feed. Figure 5-25 would result for N = 8. This also demon- 
strates that the sequential data feed is easier to implement 
in this algorithm; and, with the inclusion of the N-delay, the 
DFT algorithm of Figure 5-18 allows complete overlap of 
successive windows of data with only 7 additions required 
to set up the data rather than the 40 additions of Figure 
5-25. 
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F. CONCLUSIONS 



The convolution generator and modified convolution 

generator algorithms produce the circular convolution of 

every block of data in a continual data stream with the 

(finite) impulse response of an FIR filter. The algorithms 

are_ simpler In the number of operations to be performed and 

in topology than the definitive form of generating circular 

convolutions . 

t h 

The nr impulse response of the N outputs is a copy of 
the original FIR response circularly shifted m times and 
then delayed m sample times. Because of this the algorithm 
could be used to provide N matched filters in one algorithm, 
where the set of transmitted signals are the N circularly 
shifted versions of a sequence of length N. The maximum 
signal-to-noise ratio for a given transmitted signal will 
occur at the intersection of the diagonal line of figure 
5-15 and the output line corresponding to the circular shift 
of the transmitted signal. 

The new DFT algorithm presented here allows the program- 
ming or wiring of a DFT of length N = 2 m to be used for 
performing the DFT of lengths a lower power of two by only 
changing the parameter of length. The algorithm is designed 
for sequential data feed which allows complete flexibility 
of overlap in processing sections of data from a continual 
data stream. The zero path operations need only be performed 
when the last datum of a particular section of data has 
entered the resonator loops . No quantization errors occur 
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in the loops providing complete cancellation of previous 
data not in the section under consideration. 
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Figure 5-1 . DFT FILTER 
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Figure 5-2. Convolution of Two Sequences of Length N 
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5-3. Convolution of a Length-N Sequence with the Single Period 
Extension of another Length-N Sequence 
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Figure 5-4. Definition of Sequence (Up(n')} for N=4 
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*Im[x k (4)] T " 

Figure 5-5. State Unit Pulse Response of DFT FILTER for N=4 
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Response of DFT FILTER to 
a Periodic Input for N=4 
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Figure 5-7. Circular Convolution Generator in DFT Form 
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Figure 5-8. Circular Convolution Generator in FIR Form 




Figure 5-9. Availability Time of Convolution of Sequence 
{Up(m)} for N=8 
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Figure 5-10. Circular Convolution Generator (Definitive Form) 
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Figure 5-11. Rearrangement of Circular Convolution Generator of Figure 5-10. 
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Figure 5-13. Modified DFT FILTER 
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Figure 5-14. Modified Circular Convolution Generator 





Figure 5-15. Availability time of convolution sequence 
{Up (m) } for N=8 in modified generator 
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Figure 5-16. DFT Algorithm for N=8 
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Figure 5-17. Flow graph manipulations of DFT algorithm 
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Figure 5-18. Rearrangement of DFT Algorithm of Figure 5-16 
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Figure 5-19. Four-Delay DFT Stages; (a) with Feedback, j; (b) with Feedback, 
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Figure 5-20. Four-Delay DFT Stages with Minimum Non-Trlvial Multiplications; 
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Figure 5-21. Eight-Delay DFT Stage w:th Feedback Coefficient 
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Figure 5-22. Eight-Delay j-Feedback DFT Stage with Minimum Non-Trivial Multiplications 
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Figure 5-23. Sixteen-Delay j-Feedback CFT Stage 
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Figure 5-24. Sixteen-Delay j-Feedback i)FT Stage with Minimum Non-Trivia! Multiplications 




Figure 5-25. Parallel Input DFT Algorithm Equivalent to Figure 5-18 
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APPENDIX A 



CORRELATION IN LIMIT CYCLES 



1. Introduction 

Limit cycles occur in digital filters as a consequence 
of the non-linearity introduced by the requirement to quan- 
tize the results of multiplication. While analyzing limit 
cycles, Parker and Hess [52] showed that the quantization 
error sources exhibit limit cycles which are sinusoidal in 
nature but have a discontinuity if the principal magnitude 
is greater than the maximum quantization step size. For 
example. Figure A-la shows a limit cycle where the principal 
magnitude of the sinusoid is A while the maximum quantization 
error is M. Figure A-lb shows a limit cycle where the princi- 
pal magnitude of the sinusoid is B but B is less than M so 
no discontinuity exists. 

Returning to Figure A-la, notice that the quantization 
error departs from a pure sinusoid at most once ? always re- 
turning at the next discontinuity. Now if A had been larger 
than 3M, the next discontinuity would have to move the quan- 
tization error farther from the pure sinusoid as shown in 
Figure A-lc . In this case the next two discontinuities make 
the quantization error coincide with the original sinusoid 
again . 

From the previous discussion, define a limit cyc?.e with 
discontinuities which cause the quantization error to be at 



most once removed from the sinusoid (Figure A-la) as a case 
1 limit cycle. If the quantization error Is twice removed. 

It Is a case 2 limit cycle. A case k limit cycle Is then k 
times removed from a pure sinusoid; If no discontinuities 
exist, this is case 0. 

Now consider two quantization error limit cycles in the 
same filter. It was shown by Parker and Hess that these 
will have the same fundamental frequency but one may be 
shifted from the other by a phase 0, as x^n) and x 2 (n) in 
Figures (A-la and b). When considering two error sources at 
a time, define the limit cycle condition as case k-£ where 
one error source Is in a case k limit cycle and the other is 
In a case l limit cycle. 

It is the object of this appendix to show that quantiza- 
tion error sources are not uncorrelated, especially in the 
limit cycle condition. The correlation coefficient for limit 
cycles of the 1-0 case, and the 2-0 and 0-0 cases are con- 
sidered, assuming analog rather than digital functions for 
ease of notation. 

2. Correlation in a 1-0 Limit Cycle 

Let the limit cycle functions for a 1-0 case limit cycle 
be given by the analog equivalent of Figures A-la and A-lb 
and let 
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A sin cot 






x x (t ) 



< 



A sin cot t e {0,<J>/co} = (t x > 

or t e { ( it- <+, ) / co , (tt+<J))/(o) > - (t 3 > 

or t e {(2 tt-<J))/co,2tt/co> = {h} 

A sin cot - 2M t e { <f>/co . = { t^ } 

co 2 

^ A sin cot + 2M t e {(tt+<J>)/co, ( 2 tt_ < f> ) /to > = (t^) 



X 2 (t) = 



B sin (cot + 9) 



where 



{t i > represents the associated interval 

{-M,M> is the region of permissible values of x^(t) 

M < A < 3M 
0 < B < M 

<p = sin“ 1 (M/A) > (i.e. sin _1 (l/3) 1 <f> < tt/ 2) 
and 9 is an arbitrary phase shift. 

The correlation between these two functions is defined by 
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E[x 1 (t)x 2 (t) ] 

V2 



/ x 1 (t)x 2 (t)f T (t)dt 



°o 00 1/2 

[ / X?(t)fm(t)dt J Xp (t)f r p(t )Qt] 

— CO —00 



2tt/co 

f x-^ ( t ) x 2 ( t ) ^jfdt 



2ir/co o 2 tt/co p 1/2 

[ Q / x l (t) 5i dt Q ; x 2 (t) # dt] 
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where f,p(t) is the distribution of a random variable t, 
uniformly distributed over one period of the signals. 
Performing the integrations yields 

1 

n ( v - [t - 8r (1-r 2 ) 2 ]cos 6 

p 12^ r, °' j x 

[ir 2 - l6Trr(l-r 2 ) 2 + 8Trr 2 cos -1 ( 2r 2 -l) ] 2 
where r = M/A, i <_ r £ 1 



3. Correlation in 2-0 and 0-0 case Limit Cycles 

When the same procedure is carried out for the 2-0 case 
limit cycle, the resulting correlation coefficient is 



P 12 (r,8) 



1 1 

[tt - 8r{(l-r 2 ) 2 i- (l~9r 2 ) 2 }] 



O ^ v 



[ir 2 +l6irr 2 ( 2Tr-sin"‘ 1 'r-3sin’'' 1 '3r)-l67rr{ (1-r 2 ) 2 



+(l-9r 2 ) 2 }] 2 



For the 0-0 case, there is no dependence on the ratio r = M/A 
since then A <_ M and no discontinuities occur in x 1 (t). In 
this case the correlation coefficient is 



Pl 2 ( 8 ) = cos 0 

The correlation coefficient normalized by cos 6 is plotted 
in Figure A-2 as a function of a = 1/r = A/M for the 0-0, 
1-0, and 2-0 cases. 
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Figure A-2. Correlation in Limit Cycles 



APPENDIX B 



EXPERIMENTS IN ERROR CROSS-CORRELATION 



1. For Unquantized Signal Variable 

In order to gain an insight into the nature of the 
„ correlation between quantization error sources in digital 
filters, a simulation of the process of multiplying a single 
signal variable by two different coefficients and rounding 
the products was performed. A model of the simulation is 
shown in Figure B-l, where (u) represents a full precision, 
approximately uncorrelated Gaussian pseudo-random input with 
zero mean and unity variance. 

The object is to multiply the input signal variable by 
two coefficients, A ± and Ag, and to find the errors, e-j^ and 
e 2 , by rounding the products (indicated by the operation Q) . 
The correlation coefficient is computed by the definition 



P 12 



E[e 1 e 2 ] 




2 



(B.l) 



where E[ • ] indicates "expected value" and o q is the 
standard deviation of e^ given by 



o 



(E[e? ] ) 



1/2 



(B . 2) 



for zero mean variables. 
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To understand what occurs In the process, consider the 
error functions given in Figure B-2, where It Is assumed 
that the ratio of the coefficients Is 




and the are both positive and rounding is assumed as the 
quantization method. 

The error e^ is given by 

e i (u) = A i u - [A^]^ = A i u - k 1 (u)h (B.5) 

where h is the quantization step size and k^u) is an integer 
determined by 



k^ (u)h 





. k, (u)h 

u i-V + 

1 




(B.6a) 



which implies 




k. (u) 




(B.6b) 



or 



ki (u) 



c 



A.u 



h 



+ 



I 
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where C • ]j - greatest integer less than (•). It turns out 
then that the two functions e-^u) and e 2 (u) must coincide 
at integer multiples of p = 2h/A^ = 3h/A 2 , since e^(u) is 
periodic with period h/A 1 and e 2 (u) has period h/A 2 . These 
periods coincide when 

Jh/A x = Kh/A 2 (B. 7 a) 

(where J and K are integers) which implies that 

J/K = A^/A 2 = 2/3 (B.7b) 

by assumption. So the product e^ is periodic with period 
P. 

When e 2 (u) is plotted against e^u), the result is a 
diagram like Figures B-3(a through d) for the ratios given. 

If it is assumed that the errors are uniformly distributed^ 
between -h/2 and h/2 and that the joint probability density 
is uniformly distributed along the solid lines of Figures B-3, 
it is possible to calculate the correlation coefficient 
according to Equation (B.l). 

For example: 

Consider Figure B-3c where R is given by (B.4). The 
joint probability density is given by 



•^For a continuous Gaussian input with large enough 
variance, this assumption is well justified. 
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p(e 1 ,e 2 ) = p(e 2 ) • p(e 1 |e 2 ) 



(E.8) 



where 



P(e 2 ) 




(B. 8a) 



p (©1 e 2 ) = 



5 (e 1 -[|-(2e 2 -2h) ]) + <5 (e^C^e^h) ]) 

+ 6( e]L -[i(2e 2 )])} £ < e 2 < | 

i{<5(e 1 -[|(2e 2 -h)]) + 6 ( e n - [i( 2e 2 ) D (B.8b) 

+ 6( e;L -[i(2e 2 +h)])} < e 2 < £ 

|{6(e 1 -[i(2e 2 )]) + 6 (eyC^ 2e 2 +h) ]) 

+ 6(e r [-^(2c 2 +2h)])} < e 2 < 



and 6 ( x ) is the Dirac delta function 
Then 



E *- e l e 2 ^ 



h h 

2 2 

/ f e 1 e 2 p(e 1 ,e 2 ) de 1 de 2 
h _h 
'2 2 

h h 

? 2 

/ e 2 p(e 2 ) f e 1 p(c 1 |e 2 ) de ][ de 2 

h *1 

'2 2 
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2,2 2 e ? e 2 
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¥ 
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ITif 



(B.9a) 
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and 



h 

2 



f e 1 P(e 1 ) de i = 



h_ 

12 



(B .9b) 



so that 



12 



_1_ 

12 



(B.10) 



Bsr*f onmlng the calculation for the rest of Figures B—3 
yields for 



R = 2: p 



12 



(B.lla) 



D — 



O . 
J * 



12 



1 

3 



(B.llb) 



R = |: p 



12 



12 



(B.llc) 



R = |: P 



12 



_1 

15 



(B.lld) 



It is suggested by Equations (B.ll) that the rule for 
finding the correlation coefficient is as follows: 

(1) Write A^ and A 2 as integer multiples of their 
quantization step size, i.e. 



A ± = Jq 



(B . 12a) 



A 2 = Kq 



(B.12b) 
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Note that q is not necessarily equal to h. 

(2) Reduce J and K until they are relatively prime 
so that the ratio of coefficients becomes 



R = 




CK] 

m 



pj 

PK 



(B.13) 



where [m] means the product of the prime factors of m 
which are not prime factors of n. 

(3a) If [K]pj and [J] ^ are both odd, then 



„ _ signum [R] 

p 12 - LK JpJ LJJ p 7 



(B.l4a) 



where signum [R] is +1 for R positive, and -1 for R negative 
(R = 0 is a degenerate case). 

(3b) If [K]pj is odd and [J] ^ is even or vice 
versa, then 



„ _ 1 signum [R] 

P 12 ' ' 2 LKJ pJ UJ pK 



(B.l4b) 



The computer simulation verifies these results. There, 

and were made integer multiples of 0.01 between 0.01 

and 1.0. Ten thousand random numbers were fed into the 

simulation and the products were rounded to two decimal 

2 

places. All data for the same ratio we re averaged and 



p 

By the limits of the experiment, there was more data to 
average for low values of [J] „ and [K] , , thereby yielding 
better results. Averaging waG'nccossary due to Internal 
computer rounding ’which also occurs and interferes with 
the simulated rounding. 
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plotted In Figures B-4 through B-8. Only ratios greater 

than unity are considered and the data are divided into 

groups where the ratios are integers ([J] pK = 1), multiples 

of i ([J]p K = 2 unless K is even), etc. (Data which appears 

on a previous figure are plotted in parentheses.) 

In each figure, the values of the correlation coefficient 

p 

are also shown multiplied by ([J] v ) . When this is done 

ph 

the points fall on the two functions 1/R or -1/2R. Thus 



' l/R = CJ V [K V 



<r j V 2p = < 



or 



-1/2R = -[J] pK /2[K] pJ 



(B.15) 



and therefore 



l/[J] pK [K] pJ 



p = / or 



|^-l/2[J] pK [K] pJ 



(B .16) 



verifying Equations (B.14) except for the factor signum [R]. 
But it is obvious from Figures B-2 and B-3 that if R is 
negative, A x and A 2 are of opposite sign and the sign of the 
correlation coefficient would be the opposite of what it 
would be if A ] _ and k 2 had the same sign. (The experiment was 
performed for A-^ and A 2 both positive.) 
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2. For Quantized Signal Variable 

In Section 1, the distributions of the errors were 
assumed to be uniform. This Is because the Input had a full 
precision (essentially continuous) Gaussian distribution. 

When the input is quantized before the multiplication then 
the distributions and joint distributions become discretized 
and the results for the correlation coefficient are signifi- 
cantly different. A simulation of this case has not been 
run but the set of points in the (e l9 e 2 ) plane where the 
joint error probability density function has its support 
has been plotted for a signal quantization level of h = 1 / 16 , 
a multiplier coefficient quantization level of q = 1/8. 

These are shown in Figures B-9 . The support is indicated 
by a numeral in the graph. 

The products (EJ)h = EJ/16 and (EK)h = EK/16 are the 
errors corresponding to the multipliers = Jq = J/8 and 
A g = Kq = K/8 . If it is assumed that every point of support 
has equal probability of occurrence, then the correlation 
coefficient can be simply calculated. By summing the pro- 
ducts of EJ/16 times EK /16 for all values of EJ and EK which 
provide support and dividing by the number of points pro- 
viding support, ECe^Ogl is found. The variance of each 
error can be found by summing the squares of the values of 
EK/16 (or EJ/16) for which there is support in that row (or 
column) and dividing by the number of rows (or columns) for 
which there is support. Note that it is always true that 
if there are n (>0)' points of support in one row (or column) 
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then there are n points of support in any other row (or 
column) that contains support (as long as EK = 8 and EK = -8 
are considered the same row) . 

No attempt has been made to date to calculate the corre- 
lation coefficient for this case or to predict a formula 
for its calculation based only on knowledge of the relatively 
prime values of the multiplier coefficients, but this is 
suggested for future studies in this area. 

3. Conclusions 

a. Figures B-3 show the only possible values of and 
e 2 which can support their joint probability density function 
for the ratio given, no matter what distribution the input 
has . 

b. Equations (B.14) yield an accurate formula for the 
correlation coefficient between two sources of quantization 
error arising from the same signal variable if the distribu- 

r 

tion of that signal variable implies uniform distribution 
of the error sources . 

c. Equations (B.1*0 do not reveal any dependence of 
the correlation coefficient on the input quantization level 
or the output quantization level. It is hypothesized that 
these dependencies do exist. In addition, there may be a 
dependence of the quantization level of the coefficients and 
the common prime factors of J and K. 

d. The patterning of the distribution support suggests 
a similarity to the work of Mulcahy [55]- He examined the 
correlation between the input to a multiplier and the 
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consequent error produced by quantization. His work could 
be extended and combined with the results depicted in 
Figures (B-9) to yield insight into the nature of the 
correlation between two errors. 
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Figure B-l. Quantization error simulation model 
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Figure B-3c. Error Distribution for R 
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Figure B-3d. 
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Figure 8-4. Correlation Coefficient for Integer Ratios 
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Figure B~5. Correlation Coefficient for Ratios a Multiple 
of One-Half 



286 




\ 




K 



Figure B-6. Correlation Coefficient for Ratios a Multiple 
of One-Third 
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APPENDIX C 



VERIFICATION OF POSITIVITY OF FACTORS 
IN THE OUTPUT ERROR VARIANCE 

1 . Introduction 

In Chapter IV several expressions for the variance of 

the error at the output of a second order digital filter 

caused by quantization noise are derived. Since a variance 

is always a non-negative number, these expressions should be 

examined to ensure they yield non-negative results. 

All the expressions involve factors containing the 

coefficients of the transfer function of the digital filter as 

given in Chapter III. There are conditions on the coefficients 

in the characteristic function of the filter imposed by the 

requirement that the filter be stable. The characteristic 

function of a second-order digital filter is given by 
2 

z + az + b. For stability it is required that -1 < b < 1 
and -(1+b) < a < 1+b . There are no restrictions on the 
values which may be taken on by the coefficients c and e 
(or c’ and e’) since they do not effect the stability of 
the filter. 

2 . Positivity of the Factor, A 

The factor A, which appears in every expression for the 
output error variance is given by 

A = (1-b) [ (1+b) 2 - a 2 ] (C.l) 
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Since -1 < b < 1 it is clear that the factors (1-b) and (1+b) 
are both positive. The requirement that -(1+b) < a < (1+b) 
implies that |a[ < | 1+b | so that the factor in square brackets 
in (c.l) is positive. Therefore A is always positive inside 
the stable region. 

In addition to the positivity of A within the stable 
region of the (a,b) plane, it is obvious that A is zero on 
the boundary of that region, since the boundary is defined 
by the equations b = 1, a = 1+b and a = -(1+b). Figure C-l 
shows the regions of positivity of A. The points (a,b) = (0,-1) 
and (0,1/3) are critical points of A. The point (a,b) = (0,-1) 
is a saddle point and (0,1/3) Is a local maximum, where 
A = 32/27. At the origin, A = 1. 

3 . Positivity of the Factor, (l+b)/A 

Figure C-2 shows a plot^ of the factor (l+b)/A. It is 
clear from the previous section that this factor is always 
positive inside the stable region since A > 0 and (1+b) > 0. 
Furthermore (l+b)/A has a critical point at (a,b) = (0,0) 
which turns out to be a local minimum. The value of (1+b) /A 
at the origin is seen to be unity. 

4 . Positivity of the factor, p = [(c‘+e^) (1+b) -2ace]/A 
There are several ways to show that the factor, 

P = [(c 2 +e 2 ) (l+b)-2ace]/A > 0. First of all, note that it 



^LCDR L. Souchon Federal German Navy, 
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is only necessary to show that the numerator is non-negative 
since, in Section 2 it was shown that A > 0 inside the stable 
region. 

a. One way to show that p >_ 0 is to observe that 

(c 2 +e 2 )(l+b) > 0 (C . 2 ) 

so that it is only necessary to show that the largest value 

of 2ace is less than ( c 2 +e 2 ) ( 1+b ) . 

If the product, ce , is positive then the largest 

2ace can be is 2ce(l+b) since a < 1+b. Letting c = ye with 

2 

y > 0 forces ce = ye to be positive. Then 

(c 2 +e 2 ) (1+b) - 2ce ( 1+b ) > 0 (C.3) 

if and only if 

(y 2 +l) > 2y (C.lJ) 

and (C.4) is always true. 

For ce negative, 2ace is maximum when a = -(1+b). 

Then to show that 

( c 2 +e 2 ) ( 1+b) + 2ce ( 1+b ) ^ 0 (C.5) 

let c = -ye with y > 0. Then (C.5) is true whenever 
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(1+Y 2 ) > 2y 



as before. 

It can be noted here that the numerator of p is zero 
if c = e and a = 1+b or if c = -e and a = -(1+b) . But then 
A is zero also and p is a degenerate form. The numerator 
is also zero for any a and b inside the stable region if 
c = e = 0 . 

b. Another way to show that p > 0 is to make use of 
Section 3 by writing 

p = ( c 2 +e 2 ) [ ( 1+b ) /A] - 2ace/A . (C.6) 

For ce positive 

p > (c 2 +e 2 -2ce) (l+b)/A > c 2 +e 2 -2ce = (c-e) 2 > 0 (C.7) 

since a < (1+b) and (1+b) /A _> 1. 

For ce negative 

p >_ ( c 2 +e 2 +2ce ) ( 1+b ) /A >_ (c+e)^ 0 (C.8) 

since -a < ( 1+b ) . 

c. A third way to realize that p > 0 is to observe that 
p is a quadratic in c or e with the coefficient of the qua- 
dratic term being (l+b)/A > 1 > 0. Fixing e and finding the 
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point where the slope of p with respect to c is zero, yields 
the value of c which make p minimum for the fixed e: 



.* = 



ae 

1+b 



(C.9) 



Using c* to define the locus of points where the 

partial derivative of p with respect to c is zero, substitute 

c* into p to form p* = p| This is still a quadratic in 

I c * 

e opening upward. Taking the derivative of p* with respect 
to e then yields the value of the parameter e which yields 
the minimum of the minima 



e* = 0 



(C.10) 



which makes 




(C.ll) 



so that p is a universal minimum at c = e = 0 . Then p = 0 
unless A = 0. This case will be of no concern since when 
c = e = 0 the filter is degenerate and consists of, at most, 
a single delay element. It suffices to say that p > 0. 

5 . Nature of the factor, p+c 
The factor p+c is given by 

p+ c = ( c 2 +e 2 ) (1+b ) - 2aee + cA 12 ) 

A 
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This factor is also a quadratic in c or e opening upward so 
the same techniques as used in Section 4c can be utilized. 

Taking the partial of p+c with respect to e yields 



. * = 



ac 

1+b 



(C.13) 



for the value of e which makes p+c minimum for a fixed c 
Substituting e* into p+c to form (p+c)* = (p+c) | * and 
taking the derivative of (p+c)* with respect to c yields 



.* = 



-A (1+b) _ b^ - 1 



2[(l+b) 2 - a 2 ] 



(C.l4) 



Using this value of c in (C.13) makes the critical 
value of e 



>#* = a (k ~ 



(C. 15) 



Substituting (C.14) and (C.15) into p+c yields 



(P +c ) | c * >e ** 



_ -(1 - o ) 



(C .16) 



The minimum value of (C.l6) in the stable region is 
at b = 0 . Then 



p+c 



c=-l/2 

e=-a/2 

b=0 



= - 1/4 



(c .17) 
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Although the factor p+c is not necessarily positive 



itself, the factors in which it is used in Chapter IV are 
2(p+c) + 3 and 3p + 2c+2 and by this analysis and Section 4: 

2(p+c ) + 3 > 2.5 (C.18) 



and 



3p+2c+2 = 2 (p+c ) +p+2 > 1.5+P > 1.5 (C.19) 



since p >_ 0 . 

6 . Other Factors of Interest 

The preceding sections have been concerned with factors 
which appear in the expressions for the output error variance 
for the uncorrelated case as listed in Table 4-1. There are 

p 

factors in the general expressions for a ^ given by Equations 
(4.48) for which knowledge of positivity may be useful in 
future studies. These will be discussed here. 

a. In Equations (4.48a and c) the factor multiplying 
f^ is identical to p in Section 4 of this appendix. This 
has been shown to be non-negative. 

b. The factor multiplying f ^ in Equation (4.48a) is 

r = (c 2 b 2 (l+b)-2ce(ab 2 )+e 2 [(l+b)-a 2 (l-b) ]}/A (C.20) 
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This too is a quadratic in c or e opening upward. Taking 
the partial of r with respect to c yields 



c* = 



ae 

1+b 



( C . 21 ) 



as the value of c which minimizes r for a fixed e. Defining 
r# = r| c * and differentiating r* with respect to e yields 

e* = 0 



so that c = e = 0 yields the minimum value of r. The factor 
r is unbounded on the boundary of the stable region unless 
b=-lorc=e=0. Then it is an indeterminate form. 

c. The factor 1+a+b multiplying f^ and f^ in Equation 
(4.48d) is definitely non-negative by the requirement that 
-(1+b) £ a £ 1+b in a stable filter. 

d. The factor multiplying f^ in (4.48c) can also be 
shown to be non-negative, its minimum occuring for c = e = 0. 

e. The factor that multiplies f^p in (^»^8c) has not been 
tested completely, but it can be shown that it is a quadratic 
in c or e and that the coefficient multiplying either quadratic 
term is non-negative. It is expected that this term is 
non-negative . 

f. The factor (a) multiplying f ^ in Equation (^4.^1 8b ) 
can of course be positive or negative. 

g. The factor multiplying f ^ in (*4. 48a) can also be 
positive or negative, since the coefficient which multiplies 
the quadratic term c or e can be positive or negative. 
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Fig. C-l. Sign of A = (1-b) [ (l+b) 2 -a 2 ] in the (a,b)-plane 
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APPENDIX D 



EXPERIMENTS IN ERROR PROPAGATION 
IN SELECTED DIGITAL FILTERS 



1. Introduction 

Chapter IV presents theoretical formulae for the predic- 
tion of the variance of the output error caused by quantiza- 
tion of products in a finite precision digital filter. These 
formulae involve an infinite summation which must be checked 
for convergence. There is no doubt that the summation con- 
verges since it is assumed that the filter is stable and the 
entries of the correlation matrices are bounded. The question 
is whether in the limit this summation yields the variance 
of the output error. But a priori statistics of the errors 
are not known, particularly with respect to correlation. 

t 

In order to gain some insight into the processes involved, 
computer simulations of a TMqq filter were performed on an 
IBM 360. The objective was to test the convergence of 
Equation (^.21) as it applies to Equation (^.^8b) for various 
inputs and sets of filter coefficients. 

2. Experimental Design 

a. The Filter Simulation 

For this simulation, the canonic realization TMqq 
with d' =0 was chosen because it does not involve quantiza- 
tion errors in the output equation, it has the simplest 
expres Ion for the output error variance, and it does not 
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require any rearrangement of coefficients from those given 
in the transfer function. The state equations and the trans- 
fer function for this filter are given by: 



x-^n) 


= -ax-^Cn-l) 


+ c'u(n-l) 


(D. la) 


x 2 (n) 


= -bx^(n-l) 


+ e'u(n-l) 


(D.lb) 


v(n) = 


x-^n) 




(D.lc) 



and 



H(z) = (c> + - ? >Z -— ' ) % 1 - ( D . 2 ) 

1 + az + bz 

l 

I 

A subroutine which performed two algorithms for the 
filter was designed. One algorithm utilizes double precision 
arithmetic operations and corresponds to Figure 3-4 2, the 
other corresponds to Figure 4-la with the errors being 
generated by a quantization subroutine which rounds numbers 
to one decimal place. 

b. Statistics Collected 

The difference between the signal variables of the 
two filter algorithms yields the accumulated errors Av(n) 
and y (n) . The required statistics of these errors are 

M K 

o? ( Av) = h Z [Av(i)] 2 - [i Z Av(i)] 2 (D.3a) 

nv 11 i-1 11 i=l 
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i=l “ i=l 11 i=l 



(D.3b) 



where M is the number of iterations performed. (Writing the 
statistics as being functionally dependent on Av and y indi- 
cates the source of the statistic.) The value of M for all 
runs was 3000. 

The quanti zation errors were preserved and all auto- 
and cross-correlation functions (as defined in Appendix E) 
for these errors were calculated up to |k| = 40 i.e. 



where a, 3 = a,b,c' and/or e'. 

Furthermore, the quantization errors were summed as 
appropriate to form the composite errors e^(n) and e 2 (n) as 
shown in Figure (4-lb) and given by 



were also calculated as in (D.4) except a,B = 1 and/or 2. 




e l (n) = e (n) + e c t ( n ) 



(D.5a) 



e 2 (n) = e b (n) + e g , (n) 



(D.6b) 



The auto- and cross-correlations of these errors 



With these statistics available it is possible to 



form 
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a m . 

F(e) = R (0) + Z A 1 R (i) 
m — e ±-\ e 



m 

+ i 

1=1 



[A 1 R e (l)] t 



m=0,l, . . .40 

(D. 7) 



f* Vi 

where F^(e_) is the m “ partial sum of the right-hand side 
of Equation (4.21) and depends on the errors e. (Note that 

A 

F 0 (e) contains only R (0).) From (D.3b) it is also possible 
to form 



F(y) = A (y) - AA (y)A t (D.8) 

J *y 

which corresponds to the left-hand side of (4.21). Then 
the convergence of F m (e) to F(y) can be checked. 

2 

For the output error variance the quantities o. (y) 

uV 

2 

and cr^ v (e_) were calculated by substituting (D.8) and (D.7) 
respectively into Equation (4.48b). The convergence of 

p 2 

a^ v (e_) to a^ v (Av) can then be seen. 

The auto- and cross-correlations of the input and 
the states of the lew precision algorithm and of the state 
errors were also calculated and plotted, 
c. Filter Coefficients 

Three sets of filter coefficients were chosen to 
check the effect of changes in pole/zero locations on the 
error processes. These three sets are given in Table D-la, 
and the corresponding filters are labelled 0, 1 and 2 for 
bookkeeping purposes. Also included in Table D-l are the 
resulting parameters of the infinite precision filter for 
the given coefficients where 
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z, is the damping coefficient 

is the undamped natural frequency 

yT is the exponent of the envelope ever the 
unit pulse response 

is the frequency of the unit pulse response 
is a phasing constant given by Equations (F . 12 ) 

d. Inputs 

Three different input conditions were used to derive 
the filters, a zero mean, unity standard deviation, essen- 
tially uncorrelated Gaussian pseudo-random one, and two step 
functions, one with magnitude 1.0 and the other 0.5. The 
random input was chosen to check the shape of the state and 
state error correlations. The step inputs were chosen to 
check the response to a highly correlated input. Since this 
can also be viewed as a deterministic input, the output error 
variance should be interpreted as the mean squared error. 

- 3. Results for Random Input 

The results of the simulation for the random input are 
shown in Table D-2 where the percent normalised differences 
are shown by 

A a 2 /b = 100[a A ? v (e) - a^( Av)/a A 2 v ( Av) (D.9) 

Column (11), headed by n o^ by ( i i.66f)", indicates 
the predicted value for the output error variance given by 
Equation ( . 6 6 f ) where it is assumed that the error sources 
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are uncorrelated and uniformly distributed so that 

a 2 = h 2 /12 = (0.1) 2 /12 = 0.08333... (D.10) 

The percent normalized difference for this theoretical value 
is given in Column (12). 

All these data tend to confirm the convergence to 

the correct value. The most apparent result is that the 

rate of convergence depends very heavily on the value of 

the damping coefficient, £ . This is shovm in Column (10) 

2 

where the value of m is given for the term when a.„,(e_) 
converged to and stayed within 0.15? of the actual value. 

Thus for filter 2 the lower damping coefficient required 
that more terms be summed to get within 0.1$. It should be 
stated here that there appeared to be some divergence be- 
tween term 20 and term 40 for filter 2. This divergence 
was slow and seemed to be peaking somewhere near term 40. 

Another fact to be inferred from Table D-2 is that 
the very slight shift in zero locations between filters 0 
and 1 resulted in a detectable change in the output error 
variance. This change is more obvious for the step input 
to be discussed in the next section. This change is not 
predicted by Column (11) since the "uncorrelated" assumption 
yields formulae which do not depend on c' and e' . 

Another indication that the errors are not uncorre- 
lated is demonstrated by Figures D-2 through D-22 which show 
the correlation functions of the states and state errors. 
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(Recalling that for this filter the output is equal to state 
x 1 , figures which include this state equivalently include 
the output and figures which include the state error 
equivalently include Av.) It can be seen from the autocorre- 
lation of x 1 in Figures D-2a,~D-3a, and D-ifa that the conver- 
gence of the output autocorrelation follows the same pattern 
as the convergence of the formula for the output error vari- 
ance . That is, in Figures D- 2a and D- 3a the magnitude of the 
autocorrelation of the output enters and remains in the 
region ±0.075 at k = ±3, whereas in Figure D-^a this condi- 
tion occurs at k = ±13. The region ±0.075 is also the region 
in which non-zero values of the autocorrelation of the input 
(shown in Figure D-l) occur for k ¥■ 0. This is a measure 
of how close to uncorrelated the input is. 

These Figures also show the interesting fact that 
the accumulated state error correlations exhibit the same 
pattern as the finite precision state correlations which, 
by reference to Appendix F, can be seen to exhibit the same 
pattern as the theoretical infinite precision correlations. 
This is surprising since the quantization error sources and 
the composite errors e^ and e ^ did not exhibit any discern- 
ible pattern in their correlations (and thus were not inclu- 
ded here). There must then be some pattern in the error 
correlations which would result in the patterns of the state 
error correlations when properly combined. This combination 
might be found by examining the equation 



R (k) ^ E[y(n-k)y t (n) ] 



n-k n . . , 

= l l A n " K ” 1 E[e(i)e t (j)](A t ) n " J 

i=l j=l 

n-k n , , . 

= Z Z A n " K “ 1 R o (j-i)(A t ) n " J (D.ll) 

i=l j=l 6 

There are two more points to be made with respect 
to the results for the random input case, both involving the 
statistics of the quantization error sources. First, it 
was observed that the variance of all the error sources 

p 

was very close to h /12 as given in Equation (D.10). This 
suggests that the errors probably were distributed in accor- 
dance with the assumption of a uniformly distribution between 
-h/2 and h/2 . The mean of all the sources was zero which 
also agrees with this assumption. 

Second, when the correlation coefficient for the 
errors e c ,(n) and e e ,(n) was examined"*" it was found that it 
did not agree with that suggested by Appendix B. This con- 
firms the hypothesis of that appendix since the formulae 
given there did not hold when the input was quantized at 
about the same precision as the filter variables. 



^"Notice that the positioning of the multipliers c' and 
e' in Figure *J-1 is the same as that of the two multipliers 
in the experiment of Appendix E. 
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4 . Results for Step Input 

When the input to the filter was a step function of 
magnitude 1.0 or 0.5 some significantly different results 
were achieved. Table D-3 lists these results for mean 
squared errors rather than variances . 

Column (11), headed by 'hv^(e) by (4.91b)", indicates 
the predicted value for the mean squared output error given 
by Equation (4.91b) if the correlations are assumed to be 
constant and equal to the a posteriori value at k = 0 . 

The percent normalized differences given are calculated 
in the same way as in (D.9). 

Again the data confirm the convergence to the correct 
value. The rate of convergence also depends here on the 
damping coefficient, C, but compared to Table D-2, more 
terms are necessary because of the highly correlated nature 
of the errors. This is shown in Column (10) by the value 
of m required for convergence to within 0 . 1 %. 

Notice that the small shift in zero locations between 
filters 0 and 1 caused a significant difference in the 
resultant mean squared output error. 

The plots of the correlations for the step inputs were 
nearly constant for the states, state errors, quantization 
error sources and composite error sources. There was a 
tendency for these correlations to be less than the center 
value at |k| = 40, but this deviation from an exactly constant 
correlation is to be expected for a real situation since no 
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real power spectrum can be the ideal Dirac delta function. 

The deviation from a constant value was less than 1.0JS at 

| k | = 40. 

The mean squared error of the error sources was not, of 

p 

course, equal to h /12 for this case since this is more a 
deterministic case than a probabilistic one. Rather, it is 
most likely that the states were very nearly constant after 
the transients died out and that the error sources behaved 
similarly. Thus the only limit cycle present was a "sticking" 
effect, that is, a constant difference between the "infinite 
precision" model and the low precision model. The mean 
squared error of the error sources would then depend upon 
the value at which the states "stick". 

2 

5. Suggested Changes in the Experimental Design 

a. Filter Structure 

It is recommended that filter algorithms be developed 
for the rest of the canonic realizations of Chapter III in 
order to compare the results for different realizations. 

b. Quantization Subroutine 

Some error was introduced into the experiment by 
having a quantization subroutine which attempted to round 
to some number of decimal points while the program was run 
on a machine that operates in binary. It is a simple matter 



2 

The program used for this experiment is available from 
Dr. S.R. Parker, Electrical Engineering Department, Naval 
Postgraduate School, Monterey, California. 
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to rewrite the subroutine to round to some number of binary 
points without effecting the basic principle of the experi- 
ment. The error experienced with the present subroutine is 
probably in the fifth or sixth decimal place which can be 
significant for a case like that of filter 2 with a unity 
step input where the mean squared error had its first signi- 
ficant digit in the fifth decimal place. Additionally, sub- 
routines which perform truncation and sign-magnitude trunca- 
tion would be useful. 

c. Statistics Gathering 

The calculation of variances in this experiment was 
valid in that no variances were computed until the transients 
of the filter had decayed (after forty iterations). But 
the correlations calculated did make use of errors generated 
during the transient period. Although 3000 more iterations 
were performed, thus making the effects of the transient 
errors small, an even more accurate experiment would result 



if these errors were not included. 
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Figure D-2b. Response of Filter 0 
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Figure D-4f, Response of Filter 2 
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Figure D-4g. Response of Filter 



APPENDIX E 



SPECTRAL THEORY IN DIGITAL FILTERS 

1 . Introduction 

In fields employing analog systems, a great deal of use 
has been made of spectral theory in dealing with stochastic 
processes. In digital signal processing, although all the 
necessary tools for using spectral theory are well-defined 
[*J, 56], it has not been relied upon very heavily in this field. 
In this appendix, these tools are presented and will be 
utilized in Appendix F to gain knowledge about the spectral 
properties of the signal quantities in a digital filter. 

2 . Synopsis of Spectral Theory in the Z-Domain 
a. Power Spectrum of a Single Random Process 

Given a wide sense stationary random process, x(n), 
its autocorrelation function is defined as: 

■■ 1 

R (k) - E[x(n-k)x(n) ] k,n integers (E.l) 

Since this is a' function defined on a discrete domain 
it is possible to define its two-sided Z-transform as 

OO 

S ( z) = Z[R (k)] = Z R (k) z -k (E . 2 ) 

x x k=-~ x 

S (z) is called the power spectrum of the random process. 

X 

When normalized S (z) is called the power spectral density. 

X 
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The reason for this terminology is seen when one examines 
the relationship between R (0) and the integral of S (z) 
around the unit circle in the z-plane . 

Since R (k) is the inverse Z-transform of S (z) it 

X X 

can be written as 

R x (k) “ 2iJ / S x (z)zk ' ldz 

(All closed contour integration in this appendix is counter- 
clockwise on the unit circle.) 

Then the power in x(n) is given by 

E[x 2 (n>] = V 0) = 55j/ S x (z) ¥ 

' h A (e JuT )d U (E.H) 

— 7T 

which is the mean square value of the random process. 

Other properties of the autocorrelation functions 
are that it is an even function and that the value at the 
origin is greater than or equal to the magnitude of the 
function for any other value of the argument, i.e. 



R x (0) - l R x (k) l 



(E.5) 



Properties of the power spectrum include: 

(1) It is real and non-negative on the unit circle 
for real random processes. This follows since R (k) is real 

A 



and even. 
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(2) It Is an even function of w on the unit circle 
(z = e^ w ^). In fact. It has even symmetry with respect to 
Inversion of the argument, i.e. 

s x (z) = S x (l/z) (E.6) 

This also follows from the evenness of R (k) . 

x 

b. Cross Power Spectrum of Two Random Processes 

Given two random processes x(n) and y(n) which are 
mutually wide sense stationary, their cross-correlation 
function is given by 

i 

R v ( k ) £ E[x(n-k)y (n) ] (E.7) 

X j 

The power spectrum of this function is also given 
by its Z-transform 

00 

S (z) = E r (k)z -k (E.8) 

*y k =_a> x y 

The cross correlation function is not necessarily 
even, but it obeys the property that taking the cross corre- 
lation of the random processes in the opposite order has the 
effect of flipping the function about the vertical axis, i.e. 

V k) - v ( - k) (E - 9) 
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The magnitude properties of the cross correlation 



are : 



l R xy (k) 1 - 5 - [R x ( 0> + R y (0)] 



(E . 10 ) 



and 



|R xy OO I 2 1 R x (0)R y (0) ( E . 11) 

Property (E . 9 ) suggests that the cross power spectrum 
is not necessarily real or even on the unit circle, but that 
same property implies that 

S xy (z) * S yx (1/2> (S ' 12) 



but 



S xy (2) * S xy (1/2) 
necessarily . 

c. Power Transfer in a Digital Filter 
(1) For One Input and One Output 

Let u(n) and x(n) be two random processes which 
are related by 

oo 

x(n) = E u(n-m)h(m) (E.13) 

m=-°° 
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or 



X(z) = H(z)U(z) (E.l4) 

That is, u(n) is the input and x(n) the output 
of a digital filter with unit pulse response, h(n), and 
transfer function H(z). Then the autocorrelation of x(n) 
is given by 



R x (k) = E[x(n-k)x(n) ] 

OO CO 

= E{ £ u(n-k-m)h(m) E u(n-q)h(q)} 
m=-°° q=-°° 

CO OO 

= £ Z R (k+m-q)h(m)h(q) (E.io) 

m =_oo q=_oo u 

and the power spectrum is 



S ( z ) = E R (k) z” K 
x k=-» 

CO OO OO 

= Z Z E R (k-q+m)h(m)h(a) z - ^ 

k=_°o m=-°° q=-°° u 

OO 00 Z 

= E R u (p)z - ^ E h(m)z -m E h(q)z q 
p=_oo m=-°° q=-°° 

= S u (z)H(z)H(l/z) CE .16) 

The product H(z)H(l/z) is called the power transfer function 
for obvious reasons. 
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(2) Between Input and Output 

The cross-correlation between the input and 
output of the previous section is given by 



R (k) = 

XU 


E[x(n-k)u(n) ] 




00 

E{ E u(n-k-m)h(m)u(n)} 
m=- co 

00 

E R (k+m)h(m) (E.17) 

m=-co 



and the cross-power spectrum is 



S (z) = 
xu v 


0° 00 ^ 
E E R (k+m)h(m)z~ K 

\r=.^co rn— _oo ^ 

CO 00 

E R (p)z~ p E h(m)z m 

p=_co u m=- ro 


= 


S u (z)H(l/z) = S u (l/z)H(l/z) (e.18) 



When taken in the opposite order 



R ux (k) - 


CO 

E R (k-m)h(m) (E.19) 

m=-oo 


S ux (z) - 


S u (z)H(z) = S xu (l/z) (E . 20 ) 



(3) For One Input and Two Outputs 

Let u(n) be a random process used as the input 
to a digital filter with two outputs (or equivalently, to 
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two digital filters) and let the two outputs, x(n) and y(n) 
be related to u(n) by the transfer functions, G(z) and H(z), 
respectively. Then the cross-correlation between x(n) and 
y(n) 



R (k) = E[x(n-k)y (n) ] 

Ay 

00 oo 

= E{ l u(n-k-m)g(m) l u(n-q)h(q)} 

m=-oo q=-°° 

00 00 

= l l R (k+m-q)g(m)h(q) 

m =_oo q=_oo 

and the cross power spectrum is 



S 



xy 



(z) = 



When taken in the 

H (k) = 
yx 

S (z) = 
yx' 



00 00 00 

III R (k+m-q)g(m)h(q)z -k 

k=-~ m =_oo q=_co 

00 oo oo 

1 R (p)z” P l g(m)z m l h(q)z _q 

p =— 00 m=- co q=-°° 



S u (z)G(l/z)H(z) 



opposite order 

00 00 

1 l R (k-m+q)g(m)k(q) 

m =_oo q=_oo 

S u (z)G(z)H(l/z) = S u (l/z)G(z)H(l/z) 



= S xy (1/Z) 
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APPENDIX F 



SIGNAL CORRELATIONS FOR A TM Q q FILTER 

1 . Introduction 

The computer simulation in Appendix D utilized the TMgg 
form of a second order digital filter as described in Chapter 
III. The program produced plots of the autocorrelation and 
cross-correlation functions of various signal and error 
quantities involved in the processing of uncorrelated and 
correlated inputs. In this appendix, the spectral theory of 
Appendix E will be used to derive the expected results of 
the plots had they been performed for an infinite precision 
filter. Deviations from these predictions are discussed in 
Appendix D. 

2. Unit Pulse Responses in a TM^ q Filter 

Referring to Chapter III, the state equations for a TMqq 
filter are 



x 1 (n) 


= -ax^n-1) + x 2 (n-l) + c’u(n-l) 


(F.la) 


x 2 (n) 


= -bx^(n-l) + e’u(n-l) 


(F.lb) 


v(n) 


= x^(n) + d’u(n-l) 


(F.lc) 



Letting H^(z) and H 2 (z) be the transfer functions from 
the input to x^(n) and x 2 (n) respectively then 
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(F.2a) 



H 1 (z) 

H 2 (z) 



[c' + e'z ^ ]z ^ 

IT 1? 

1 + az + bz 

[e' + (ae '-bc ’ ) z -1 ]z 

-1 -2 

1 + az + bz 



-1 



(F.2b) 



Note that for d' = 0, H^(z) = H(z), the transfer function 
of the filter, but for d' = 1, H^(z) is as given in (F.2a) 
but the transfer function to the output is 



H(z) = z" 1 + I^- + e 



1 + az 



Z-^Z' 1 -1 , TI („\ 

- _2 =2 + H ^ z ) 



(F.3) 



+ bz 



Now consider the Z-transforms of the two functions 



g 1 (n) = 



e~ yn ^cos gnT 



0 



n > 0 



n < 0 



(F.iJa) 



e~ yn ^sin gnT 



g 2 (n) = 



n > 0 
n < 0 



(F.Hb) 



which are 



p ,„s. _ 1 - e yT cos BTz - " 1 

1 1 - 2e“ yT cos BTz -1 + e” 2YT z“ 2 



(F.5a) 



r /„\ _ e yT sin BTz 

u 2' ,z; -vT II -2 yT -2 

1 - 2e Y cos BTz + e r z 



(F.5b) 
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Any transfer function of the form 



(c n +c 0 z" 1 )z" 1 

G(z) = 

1 + az“ + bz 



can be written as 



(F.6) 



G(z) = c 1 CG 1 (z) + kG 2 (z)]z 1 



(F.7) 



by writing 

, , ~lv -1 

(c,+c 0 z )z 

G( z ) = 

1 + az + bz 

c 1 (l - e~ Y ^cos gTz - ^ + ke -Y ^sin 6Tz~^)z~‘ 
1 - 2e“^cos gTz -1 + e - ^ YT z"" 

Equating coefficients it is found that 



c (1 + 3 Z -1 )Z -1 
1 C 1 
1 + az -i + bz 



= c i [G^ ( z)+kG 2 ( z) ] 



Y T = - | In b 



(F.9a) 



6T = cos -1 ( -a/2~\/~b" ) 

— + e" YT cos 6T 

k = 

e YT sin 6T 



2c 2 - ac^ 



2c, (b - a 2 ) 1//2 Vb~ 
1 T 



(F.9b) 



(F.9c) 



Here 6 is the discrete response frequency and y = £u n 
where £ is the damping coefficient and is the undamped 
natural frequency of the filter. 
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Taking the Inverse Z-transform of G(z) then 



g(n) = c 1 [g 1 (n-l) + kg 2 (n-l)] 

- c^e“"^ n “‘^ [cos 6(n-l)T + ksin 8(n-l)T] n >_ 1 

0 n £ 0 

(F . 10) 



For the transfer functions of Equations (F.2), the unit 
pulse response of the states of the TM Q q filter are 



h^(n) = 



c , e -Y(n-l)T[ cos g( n _ 1 ) T + kiS i n $( n -l)T] n > 1 



n ^ 0 
(F.lla) 



h 2 (n) = 



i e " Y (n " 1)T [ C os 8(n-l )T + k 2 sin B(n-1)T] n > 1 



n <_ 0 
(F.llb) 



where 



k„ = 



2e 1 - ac 



2c* (b - a 2 ) 1/2 Vl 

T 



(F . 12a) 



k~ = 



ae 1 - 2bc 



2e 1 (b - a 2 ) 1/2 VT 
~T 



( F . 12b ) 



and y and B are given by (F.9a and b). 
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3. State Response to an Uncorrelated Input 



When the input, u(n), to a digital filter is a wide 
sense stationary uncorrelated random process its autocorrela- 
tion function is given by 



R u (k) = E[u(n-k)u(n)] = 



where Nq is the mean square power in u(n) . Its power 
spectrum is then 



N, 



k = 0 



k ? 0 



(F.13) 



S u (z) = 



Z[R u (k)] = 



l 

k=_oo 



R u (k)z 



-k 



= N, 



(F.14) 



Referring to Appendix E, the power spectra for x., (n) 
and x 2 (n), the states of the TM^ filter, are 



S (z) = S u (z)H 1 (z)H.(l/z) = NqH.CzJH.CI/z) 

S x Xz (z) = S u (z)H 1 (1/z)H 2 (z) = N q H 1 (1/z)H 2 (z) 



(P.15) 

(F.l6) 



S (z) = S (z)H.(z) = N n H, (z) (F .17) 

ux^ u i 0 i 

For d’ = 0, the output power spectrum is given by (F.15) 
with i = 1 since then v(n) = x^(n). This is also the cross 
spectrum between the output and x-^(n) . The cross spectrum 
between the input and output is (F.17) with i ~ 1 and the 
cross spectrum between the output and x 2 (n) is given by (F.l6) 



381 



The auto- and cross-correlations involved are by defini- 
tion the inverse Z-transform of the spectra. But observe 
that each of the spectra is the product of two functions in 
the z-domain. By the convolution theorem this corresponds 
to the convolution of the inverse Z-transforms of the two 
functions being multiplied. If the inverse Z-transform of 
H^(z) is hm(k), then the inverse Z-transform of H^(l/z) is 
h i (-k) so that 



(F.ll) require that h^(k) be zero for k <_ 0, (F.18) can be 
written as 




oo 



CO 



E h. (m)h . (-(k-m) ) = E h. (m)h, (m-k) 



(F.18) 



which is the correlation of h i (k) with h^ (k) . Since Equations 



OO 




E h . (m)h . (m-k) 
m=l 1 J 



k < 0 



OO 



_ m=k+l 



E h. (m)h.(m-k) 
=i/+i 1 J 



k > 0 



OO 



E h. (m+l)h . (m-k+1) 
m=0 1 J 



k < 0 



OO 



E h (m+l)h . (m-k+1) 
m=k 1 J 



k > 0 (F.19) 
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Considering the functions of Equations (F.ll), then 



N 0 h ± (m+l)h .(m-k+1) = K ± 2n " k ^ T [cos 0mT + k ± sin 0mT] 



[cos S(m-k)T + kjSin B(m-k)T] 



(F . 20) 



since the indices m+1 and m-k+1 are at least unity for the 

values of m and k in (F.19), and is a constant involving 

Nq and the scaling constants c' and e' in (F.ll). 

By a change of variable in the second summation of (F.19) 

the correlation function R. . (k) becomes 

ij 



R ij' k) -J 



K. e” y l k l T S e ° T-T [ccs BmT + k.sin Bmt] k < 0 

1 . J r\ — 



m=0 



[cos 0(m-k)T + k^sin B(m-k)T] 



I I 00 

K. . e " y l k l T S e" 2YmT [cos 0(m+k)T + k.sin 0(n+k)T] 
1J m=0 1 



[cos 0mT + k.sin BmT] 
J 



k > 0 
(F. 21) 



It can be shown by the use of trigonometric identities that 



R . . ( k ) = K. .e" Y l k l T {<t ii cos 0kT + ijn . sin B | k | T } (F.22) 
i J 1 J v xj 
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where 



U 



, (1-k k ) (1-e 2YT cos26T) + (k.+k,)e“^’ T:: sin2BT 1+k.k. 

k LJ - — kJ- «• 



l-2e“ 2YT cos26T+e ^ 



1-e 



-2yT 



(P.23) 



* 



ij 



, (k +k ) ( 1-e 2YT cos2BT)-(l-k,k.)e 2YT sin2BT k,-k. 

= — r i ^ r _± H -J — i 

2 l-ae-^cosZST+e -^ 1 'l- e - 2r; ' 



(F . 2 1 ! ) 



using the plus sign for k ^ 0 and the minus for k _< 0 . Notice 
that for i = j, R^.Ck) is an even function of k and its maxi- 
mum occurs at k = 0, so that R^Ck) is a valid auto-correlation 
function . 

The point to be made here is that R^.Ck) is bounded in 
magnitude by k^ ^ (4>^j + ^ij) e - ' 1 ^" 1 . For i = j, this bound 
is unambiguous. But for i / j the choice of the sign in the 
expression for must be the one which maximizes . None- 
theless, this bound has the same form for every correlation 
function involving the states and output of the infinite 
precision TM^q filter with d* = 0 when the input is an 
uncorrelated wide sense stationary random process . A similar 
analysis would show that this is true for any filter realiza- 
tion with the same characteristic function. The only excep- 
tion is that at the origin, R y (k) , the output autocorrelation 
will have an additional constant term when d (or d’) =1. 

This occurs because d = 1 implies a transfer’ function of 
unity between the input and output in parallel with the 



rest of the transfer function. Then the input autocorrela- 
tion pulse at k = 0 is transferred directly to the output. 

To show this, let 

= H(z) = 1 + G( z) (F.25) 

where G(z) is in the form of (F.6). Then 

H(z)H(l/z) = 1 + G(z) + G(l/z) + G(z)C-(l/z) (F.26) 



and 



N 

Hy(k) = H(z)H(l/z)dz/z 



for an uncorrelated input so that 



(F.27) 



1 + 2g( 0 ) + Z g^(m) 
m=0 



irV k) - J 

o 



g(-k) + Z g(m) g(m-k) 
m=0 



k = 0 

k < 0 (F.2S) 



g(k) + Z g(m) g(m-k) k > 0 

m=k 



since g(n) = 0 for n < 0. So for k / 0 R (k) is bounded 

I k I T 

by a constant times e 1 1 but at the origin there is an 

additional [1 + g(0)]N term. 
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. State Response to a Highly Correlated Input 



In the previous section an uncorrelated input was used 
to determine the correlation functions of the states of the 
filter. The input power spectrum in that case was found to 
be a constant. At the other extreme, suppose the auto-corre- 
lation of the input is a constant. This is an extreme case 
because as noted in Appendix E the magnitude of the auto- 
correlation function at the origin is greater than or equal 
to the magnitude at any other value of the argument. It 
will be shown here that this kind of input yields correlation 
functions of the states which are also constant. 

a. The Power Spectrum for a Constant Autocorrelation 
The first problem involved in finding the power 
spectrum of a constant autocorrelation function is that a 
constant discrete function is not absolutely summable so that 
the application of the two-sided Z-transform may not be 
valid. This difficulty can be removed by letting the 
constant autocorrelation be given by"*" 

R (k) = lim e*" Y ^ a >_ 0 (P.29) 

u a+0 

It Is now possible to apply the two-sided Z-transform 
to 



■'"This method was suggested by Prof. Craig Comstock, Cdr. 
USNR, of the Mathematics Department, Naval Postgraduate School. 
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(P.30) 



R“(k) * e-“l k l 



since this function is absolutely summable . 
The power spectrum of R^(k) is 



S,“(z) = Z[R®(k)] = Z e " a l k l z “ k 

k=-°° 



u 



u 



OO 00 

- -ak -k . -a z v -ak k 

Z e z + e E e z 



k=0 



k=0 



. -a 

1 + e z 



, -a -1 , -a 

1 - e z 1 - e z 



(F. 3D 



as long as z is in the annulus e“ a < |z| < e a . Then by 
collecting terms 



, -a ou 
z(e - e ) 

\x t -a x , ot\ 

(z - e ) ( z - e ) 



S, a (z) = 



(F.32) 



Letting S^(z) be the analytic extension of the power spectrum 
and applying the inverse Z-transform, then by residue theory 



i“(k) = ^ i s“(z)z k 32 = 0 ■ (e ' C ‘ 

u v ' 2 ttj u z 2ttj -a,. 



R 



2ttj u 



(z - e u ) ( z - e a ) 



dz 



-ak 



ak 



k > 0 



k < 0 



= e 



-a | k | 



(F . 33) 
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so that this is a consistent transform pair. Notice that in 



the limit the power spectrum S^(z) becomes the Dirac delta 
function at z = 1, since 



S (z) = Z[R (k)] = Z [lim R^OO] 

a+0 u 

00 ii 00 i i 

— v -| -? -a k ~k *, -a k -k 

= l lim e 1 *z = lim Z e 1 'z 

k=0 a+0 • a+0 k=0 

00 z = 1 

= lim s“(z) = = 6 ( z— 1 ) (P.34) 

a->-0 0 otherwise 

This yields the transform pair 



R (k) S (z) 

u u 

1 <5(z-l) 



(F • 35a) 



1 

0 



k = 0 
k ? 0 



1 



(F. 35b) 



b. State Correlations 

It is now possible to find the correlations among 
the states and output for a constant input autocorrelation. 
In any second order digital filter, there are two transfer 
functions between the input and the states and one between 
input and output. Let them be denoted by H^(z), H^Cz) and 
H^Cz) respectively. 

From Appendix E, the power spectra for the states 
and output are given by 
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(F. 36) 



S iJ (z) = S u (z)H 1 (l/z)H (z) 

where i,j = l, 2 ,v and S^Cz) is the power spectrum of the 
input. For the case of the highly correlated input, let 

S, ,(z) = lira S, a .( z) = H. ( l/z)H . ( z) lim S a (z) CF. 37) 

1J a -»-0 J i j a ^. 0 u 

ct 

and consider S^.(z). The corresponding correlation function 
is then 

R ij (k) = 2 iJ $ H i (l/z)H j .(z)S^(z)z k ^ CF. 38) 

Now for any second order filter the product H.(l/'z)H.(z) 

J 

can be written as 

H i (l/z)H j .(z) = K ij (z) (F. 39 ) 

(z-z 1 ) (z-z 2 ) ( z-z” )(z-z“^) 

where z, and z 0 are the poles of the filter and K 4 .(z) is 

Id JL J 

a polynomial in z. Then the integrand in (F. 38 ) becomes 

I(z) = H. (z)H.(z)S a (z)z k_1 
1 J U 

.. / \ / ™ct a N k 

K ij . (z) (e - e )z 

( z-z 1 ) (z-z 2 ) ( z-z” 1 ) ( z-z^ 1 ) ( z-e” a ) ( z-e a ) 

(F.iJO) 
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For k >_ 0, it is possible to write the partial expan- 



sion of (F.iJO) as 

A (k) B (k) B (k) C n (k) C p (k) 

_± + -k — + Jl — ^ + -£ ] 

z-z~ -1 -1 -a a 

2 z-z^ z-z 2 z-e z-e 

(F.i]l) 

where C 1 (k) and C 2 (k) contain the factor ( e~ a -e a ) in their 
denominators while the other residues do not. Letting 

C 1 (k) = C^(k)/(e" a -e a ) (F.H2) 

then by residue theory 



Rjj ( k ) = [A 1 (k)+A 2 (k) ](e a -e a ) + C^(k) k > 0 

(F.H3) 



I(z) = (e~ a -e a )[ 



A n (k) 



Z-Z- 



so that 



R. . (k) = lim R. a . (k) = lim C’ (k) = lim I(z) (z-e~ a ) I 

a+0 a-*-0 1 cx+0 lz=e 



= lim e ak H. (eo)H.(e a ) 

a+0 1 J 



= H jL (l)Hj(l) k > 0 (F.M) 
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It can be shown that Equation ( F . 4 4 ) also holds for k < 0 
except that the last limit is 

lim {e” ak H (e a )H (e _a ) + (e ak - e" ak )H, («)H (0)} 
a-»-0 x j i J 

= H i (l)H J (l) (P.H5) 

as long as the filter has no pole at z = 0. 

The final result is that if the input to a digital 
filter has a constant autocorrelation function equal to N^ , 
then the auto- and cross-correlations of the states and output 
is given by 

R ±J (k) = N^CDH (1) (F.46) 

(where i,j = l,2,v), which is also a constant correlation 
function. It can also be shown that the cross-correlation 
between the input and one of the states or the output is 

00 

R ,00 = N n H,(l) = N n I h. (n) (F.H7) 

U1 u x u n=0 

Notice that the cross- correlation functions are even func- 
tions in this case. Notice also that the quantity (1) is 
the sum (discrete integral) of the unit pulse response, h^(n), 
of the associated state or output . 
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APPENDIX G 



TRANSFORMING (m) th ORDER COMPLEX FILTERS 
INTO (2m) th ORDER REAL FILTERS 



In Chapter V, the DFT Filter is an algorithm which con- 
tains first order complex coefficient digital filter sections. 
These filter sections are purported to be digital resonators 
at certain frequencies around the unit circle in the Z-plane . 
But a resonator is expected to be a second order real coeffi- 
cient algorithm. It is shown here that the equivalence is 
justified. 

First consider the first order digital filter using 
complex arithmetic with one input and one output (see Figure 
G-l) : 



x(n) = ax(n-l) + bu(n) 



(G.la) 



y(n) = cx(n) + du(n) 



( G . lb ) 



where a, b, c, d, u, x and y are complex numbers. Let 



x(n) = x^n) + jx 2 (n) 



(G.2a) 



y (n) = y-^n) + jy 2 (n) 



(G.2b) 



u(n) = u-^n) + ju 2 (n) 



(G.2c) 
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a = 8l ± + Ja 2 (G. 2d) 
b = bi + jb 2 ( G . 2e ) 
c = c 1 + jc 2 (G.2f) 
d = d x + jd 2 (G . 2g) 



where subscript (1) represents the real part of the complex 
number and subscript (2) is the imaginary part. Then Equa- 
tions (G.l) may be written 



x x (n) + jx 2 (n) = a 1 x 1 (n-l) - a 2 x 2 (n-l) + j[a 1 x 2 (n-l) 



+ a 2 x 1 (n-l)] 



+ b 1 u 1 (n) - b 2 u 2 (n) + jCb^Cn) + b 2 u 1 (n)] 



(G.3a) 

y-^n) + jy 2 (n) = c^Cn) - c 2 x 2 (n) + j[c 1 x 2 (n) + c^Cn)] 

+ d 1 u 1 (n) - d 2 u 2 (n) + jtd^^n) + d 2 u 1 (n)] 

(G. 3b) 

Equating the real and imaginary parts of Equations (G.3) and 
using state variable form yields: 
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’x 1 (n)' 




x 2 (n) 





a i " a 2 



x 1 (n-l) 



l_x 2 (n-l) 



b l " b 2 



L b 2 





u 1 (n)' 




u 2 (n) 



(G.4 a) 



'y 1 (n)‘ 

. y 2 (n) - 




"x 1 (n)“ 




[ 

C\J 

\ 

i — 1 

V 




~u 1 (n)‘ 




+ 








1 

X 

rv> 

9 s 

1 




I 

< — 1 

OJ 

i 




u 2 (n) 








(G.4b) 



These can be recognized as the equations for a second order 
digital filter using scalar arithmetic with two inputs and 
two outputs (see Figure G-2). The characteristic equation 
of this second order system is: 



z 2 - 2a 1 z + a 2 + a 2 = z 2 - 2R e Ca]z + |a| 2 (G.5) 



The pole locations are therefore z = a and z = a* where * 
indicates complex conjugate. Note that if (a) is a real 
number then a 2 = 0 and a = a* and this second order filter 
is actually two uncoupled first order filters each with a 
pole at z = a = a^. 

The general result is that an m bb order digital filter 
using complex arithmetic with p complex inputs and q complex 
outputs can be transformed into a (2m) bb order digital filter 
using scalar arithmetic with 2p scalar inputs and 2q scalar 
outputs. This is shown in the following development. 
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t h 

For the in order complex filter we have the equations: 



x(n) = Ax(n-l) + Bu(n) 


(G.5a) 


y(n) = Cx(n) + Du(n) 


(G.5b) 



where 



x(n) 


is 


an m x 1 column 


vector 


u(n) 


is 


a p x 1 column 


vector 


y(n) 


is 


a q x 1 column 


vector 


A is 


an 


m x m matrix 




B is 


an 


m x p matrix 




C is 


a 


q x m matrix 




D is 


a 


q x p matrix 





e .g. 





1 

"c 

1 — 1 
X 

i 




a n a l 2 • • • a im 


x(n ) = 


x 2 (n) 

• 

• 


A = 


a 21 a 22 

• • • 

• • • 




x (n) 

L m J 




• • • 

a • ■ • • a 

ml mm 



( G . 6 ) 



Let 

Rx^Cn) be the real part of x^(n) 
Ix,(n) be the imaginary part of x 1 (n) 
Ra^ be the real part of a^. 
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Ia. . be the imaginary part of a. 
etc . 



Let RA = (Ra ±j .) IA = (Ia ±J ) Rx = (Rx ± ) 
i.e. RA is the matrix whose elements are Ra^ . 
real and imaginary parts in (G.5) and forming 
vectors and matrices as follows: 



etc . 

Then equating 
partitioned 



Rx(n) 




~RA 


: -ia' 

i 

i 


Ix(n) 




i 

IA j RA 


" Ry(n)~ 




"rc 


-IC ' 


Iy(n)_ 




IC 


RC 



Rx(n-I) 




” RB J 

i 

i 


-IB 




Ru(n) 


Ix(n-1 ) 


1 


IB j 


RB 




Iu(n) _ 



(G.7a) 



Rx(n) 




RD ! -ID 
1 
1 


Ru(n) 


Ix(n) 




1 

ID ! RD 


Iu(n) 



(G.7b) 



Another form of these equations is found by replacing each 
complex element in each complex vector and matrix with its 
equivalent scalar vector or scalar matrix defined as 



x ± (n) 



a 






Rx ± (n) 








4 *i 


(n) 


Ix i (n) 






Ra « 


-Ia. . 1 
ij 








£ 


Ia « 


Ra ij. 






(G . 8a) 



(G. 8b) 



etc . 
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Then Equations (G.5) become 



X-^Cn) 



X 2 (n) 



^11 






a 2 i 



il2 



A, 

— lm 



—22 

X 

. I 






ini 1 



I A 
-mm 



X-^Cn-l) 



X 2 (n-1) 



j X (n-i) 

L“ m . 



- 11 ' &12 
J 



— 21 ^ -22 



B 



!p 



ini' 



'B 

! — nrp 



U^n) 



(G . 9a) 



Yi(n) 



-•li ! —i2 i 



i n 
i — lm 



— 2 ! ! 



r 



iq(n) 



C 

-ql 



qm 



1 


X x (n) 






• 






• 

• 

X (n) 


+ 




— m 





En 



E21 



Eia! '• 



I 

*V 

ql! 



D- •. i 



IP 



• 1 D 

I qp 



U (n) 
-P 



(G.9b) 



Equations G.7a and G.9a each represent 2m real state equations 
with 2p real inputs. Equations G.7b and G.9b each represent 
2q real output equations with 2p real inputs. 
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Figure G-2. Second Order Real Coefficient Digital Filter 
Equivalent to Figure 6-1. 
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APPENDIX H 



COMPLEX INTEGRATION INVOLVING A 
COMPLEX DIRAC DELTA FUNCTION 

From the results of Appendix F it is possible to propose 
the following theorem. 

THEOREM : If H(z) is an analytic function on the unit circle, 

then if H(z) has only simple poles inside the unit 
circle 

6 H(z) 6 ( z-1) ^ = 2irjH(l) (H.l) 

z 



where the integral is taken counter-clockwise around the unit 
circle and <5 (z-1) is a Dirac delta function defined by 



6 (z-1) - lim 
ct+0 + 



z(e~ a - e a ) 
(z-e" a )(z-e a ) 



“ z = 1 

0 z 1 1 



(H . 2 ) 



PROOF: Consider the integral 



1(a) 



H( z) (e~ a -e a ) 

( z-e*" a ) ( z-e a ) 



dz 



a > 0 (H.3) 



If H(z) has k simple poles inside the unit circle, then by 
residue theory 



1(a) = 2tt j (e~ a -e a ) E R, + 2irjH(e" a ) (H.l|) 

i=l x 
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where R. is the residue of at the i*'* 1 pole 

1 ( z-e~ a ) ( z-e a ) 

of H(z) inside the unit circle. Notice that does not 

have the factor (e"" a +e a ) in its denominator unless H(z) 

has a pole at e“ a . To avoid any conflict assume that the 

pole with the largest modulus inside the unit circle has 

|z| = b, and restrict a such that 0 < a < -In b. Under 

this condition 



lim 1(a) = 2 ttjH(1) (H.5) 

a-*0 



On the other hand, by (H.3) 



lim 1(a) = 6 H(z)6(z-1) — (H.6) 

a-0 z 



QED 

The previous theorem leads to the hypothesis that if 
ip is on the closed contour C, and H(z) is analytic on C and has 
only simple poles inside C, then 

c 6 H( z) 6 ( z-ip) — = 2irjH(^) (H.7) 

z 

where the integral is taken counter-clockwise around C. 
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